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ABSTRACT 

By means of two dimensional, general-relativistic, magneto-hydrodynamical simulations we 
investigate the oscillations of magnetized neutron star models (magnetars) for one particu- 
lar dipolar magnetic field configuration including the description of an extended solid crust. 
The aim of this study is to understand the origin of the quasi-periodic oscillations (QPOs) 
observed in the giant flares of soft gamma-ray repeaters (SGRs). We confirm our previous 
findings which showed the existence of three different regimes in the evolution depending 
on the magnetic field strength: (a) a weak magnetic field regime B < 5 x 10 13 G, where 
crustal shear modes dominate the evolution; (b) a regime of intermediate magnetic fields 
5xl0 13 G< B < 10 15 G, where Alfven QPOs are mainly confined to the core of the neutron 
star and the crustal shear modes are damped very efficiently; and (c) a strong field regime 
B > 10 15 G, where magneto-elastic oscillations reach the surface and approach the behav- 
ior of purely Alfven QPOs. When the Alfven QPOs are confined to the core of the neutron 
star, we find qualitatively similar QPOs as in the absence of a crust. The lower QPOs asso- 
ciated with the closed field lines of the magnetic field configuration are reproduced as in our 
previous simulations without crust, while the upper QPOs connected to the open field lines 
are displaced from the polar axis. The position of these upper QPOs strongly depends on the 
magnetic field strength. Additionally, we observe a family of edge QPOs and one new upper 
QPO, which was not previously found in the absence of a crust. We extend our semi-analytic 
model to obtain estimates for the continuum of the Alfven oscillations. Our results do not 
leave much room for a crustal-mode interpretation of observed QPOs in SGR giant flares, 
but can accommodate an interpretation of these observations as originating from Alfven-like, 
global, turning-point QPOs (which can reach the surface of the star) in models with mean 
surface magnetic field strengths in the narrow range of 3.8 x 10 15 G < B < 1.1 x 10 16 G (for 
a sample of two stiff EoS and various masses). This range is somewhat larger than estimates 
for magnetic field strengths in known magnetars. The discrepancy may be resolved in models 
including a more complicated magnetic field structure or with models taking superfluidity of 
the neutrons and superconductivity of the protons in the core into account. 



1 INTRODUCTION 

The detection of quasi-periodic oscillations (QPOs) in the giant 
flares of Soft Gamma-ray Repeaters (SGRs) has raised hopes of 
drawing conclusions about their interior structure, and thus to in- 
crease the understanding of the equation of state (EoS) of matter 
at supranuclear densities. SGRs are a class of magnetars, i.e. very 
compact objects with very strong magnetic fields ( [Duncan &| 
[Thom pson 1992 ). They show repeated activity in the soft gamma- 
ray / hard X-ray spectrum, and as of today giant flares with energies 
between 10 44 — 10 46 erg/s have been detected in three of these 
objects. In the decaying X-ray tail of two of these events a num- 
ber of QPOs have been observed (see |Israel et aL 2005 ; W atts &| 
Strohmayer 2007 for recent reviews). The frequencies of the QPOs 
are roughly 18, 26, 30, 92, 150, 625, and 1840 Hz for the outburst 



of SGR 1806-20 and 28, 53, 84, and 155 Hz for SGR 1900+14. Re- 
cently, El-Mezeini & Ibrahim (2010 ) re-analyzed the Rossi X-ray 
Timing Explorer (RXTE) observations and claimed that they dis- 
covered 84, 103, and 648 Hz QPOs in the normal bursts of SGR 
1806-20. With a different method [Hambaryan et al] ( [20TT] ) have 
also found new QPOs in the data of the SGR 1806-20 giant burst 
at frequencies 16.9, 21.4, 36.4, 59.0, 116.3Hz. If confirmed, these 
are truly fascinating discoveries as the methods employed will al- 
low one to increase the number of extracted QPO frequencies. 

The interpretation of the observed QPOs in terms of oscilla- 
tions of the magnetar itself seems very promising. It may allow for 
insight into the properties of these objects and constrain the EoS 
above nuclear matter density. The first step towards such magnetar 
asteroseismology would be the identification of the modes of the 
star which have frequencies in the right range and which could be 
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observable during an outburst. For some time after the discovery 
of the QPOs in SGR 1900-16 in 1998 the main focus was on the 
torsional shear oscillations of the solid crust because their frequen- 
cies in the range between 10s of Hz for nodeless modes and kHz 
for n = 1 modes match the observed frequencies cited above (see 
| Duncan| 19981 |Strohmayer & Watts|2005l|Piro|2005[|Sotani et al.| 
|20071|Samuelsson & Andersson|20071 |Steiner & Watts|2009| and 
references therein). From energetic considerations it is very likely 
that these oscillations are excited during a SGR outburst ( |Duncan| 
1998 Levin & van Hoven 201 1). Moreover, torsional shear oscilla- 
tions couple preferably to the exterior magnetosphere jBlaes et aT] 
1989), where the emission in form of a trapped fireball is supposed 
to occur ([Thompson & Duncan 2001). Therefore, there is a natural 
channel of how these oscillations may influence the X-ray signal 
emitted during the flare. However, despite some recent improve- 
ments of the models (Steiner & Watts 2009), the order of the fre- 
quencies of successive shear modes does not allow for a complete 
interpretation of all observed QPOs. 

Since the publications by Levin (2006) and Glampedakis & 
Andersson ( 2006) another possibility has been investigated. These 
authors show that the crust-core coupling due to the extremely 
strong magnetic field present in magnetars may have significant 
impact on the shear oscillations of the crust. It was shown in simpli- 
fied toy models that the shear modes can be damped very efficiently 
into a magneto-hydrodynamic (MHD) continuum of Alfven oscil- 
lations existing in the core of the neutron star. This idea stimulated 
further interest in the direction of magneto-elastic oscillations by 
more detailed toy models ( Levin 2007 ; Lee 2007 2008). In partic- 
ular,[Levin ( 2007) showed that due to the coupling through the crust 
long-lived QPOs can be produced at the turning points or edges of 
the continuum. 

Before studying magneto-elastic oscillations in realistic sce- 
narios it is necessary to understand purely Alfven oscillations 
of neutron stars in General Relativity. Studies using general- 
relativistic MHD models without taking an extended crust into ac- 
count revealed two families of long-lived QPOs in the continuum 
formed by Alfven oscillations of dipolar magnetic field config- 
urations jSotani et al1|008b| |Cerda-Duran et aL|[2009| |Colaiuda| 
|et al.||2009) . [Cerda-Dur an et al.| ( |2009| > derived a semi-analytic 
model based on standing waves in the short-wavelength approx- 
imation, which reproduces the MHD numerical results with very 
good agreement. The first family of QPOs, termed upper QPOs, is 
related to the turning point of the continuum at the open field lines 
close to the polar axis (see Fig.[TJ. The second family, the lower 
QPOs, can be found at the turning point in a region of closed field 
lines near the equator. In this model the upper QPOs have their 
maximum amplitudes at the surface of the star and are therefore 
candidates to explain the observed QPOs. The different overtones 
have frequency ratios given by integer numbers, making them very 
attractive to explain some of the observed frequencies, such as the 
30, 92 and 150 Hz QPOs in SGR 1806-20. 

Very recently a few groups have published results consider- 
ing magneto-elastic oscillations in magnetar models, pan Hoven"&| 
Levin ( 2011) have presented a non-relativistic model and shown 
that the oscillatory spectrum can be influenced significantly by the 
physical properties of the model. For example the presence of an 
entangled magnetic field allows for the existence of discrete Alfven 
modes, or the decoupling of neutrons in the core from Alfven waves 
can change the frequencies of the continuum. In a different ap- 
proach, based on the coupling between a number of linear one di- 
mensional wave equations in the core and a two dimensional wave 
equation in the crust, Colaiuda & Kokkotas ( 2011) suggested the 
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Figure 1. Upper panel: Illustration of the dipolar magnetic field configura- 
tion including field lines which close inside the neutron star (red lines) and 
open lines which extend to the exterior (blue and black lines). Lower panel: 
The frequencies of the corresponding field lines. The ensemble of the fre- 
quencies of all field lines forms a continuum with edges and turning points 
as indicated with the arrows. We call the QPOs related to the turning point 
of the open field lines upper QPOs, the QPOs related to the turning point of 
the closed field lines lower QPOs, and the QPOs related to the edge of the 
continuum at the last open field line edge QPOs. 



existence of global, discrete Alfven modes in the gaps between con- 
tinuous parts of the spectrum, as was also suggested in | van H oven 
|& Levin| ( [20TT] ).Inthat case, the shear modes would not be damped 
resonantly and they could be observed as QPOs. These oscillations 
are called gap modes. 

In a previous paper ( Gabler et al. 2011), we presented a re- 
alistic model of magneto-elastic oscillations in magnetars, finding 
that crustal shear oscillations, often invoked as an explanation of 
the QPOs seen after giant flares in SGRs, are damped by resonant 
absorption on timescales of at most 0.2s, for a lower limit on the 
dipole magnetic field strength of 5 x 10 13 G. At higher magnetic 
field strengths (typical in magnetars) the damping timescale is even 
shorter. Our findings thus exclude torsional shear oscillations of the 
crust from explaining the observed low-frequency QPOs. In addi- 
tion, we found that the Alfven QPO model is a viable explanation 
of observed QPOs, if the dipole magnetic field strength exceeds 
a minimum strength of about 10 15 G. Then, predominantly Alfven 
QPOs are no longer confined to the fluid core, but completely domi- 
nate over shear oscillations in the crust region and have a maximum 
amplitude at the surface of the star. 

In the following, we present the details of the extension of a 
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previous model by Cerda-Duran et al. (2009), where we performed 
two dimensional MHD simulations of general-relativistic magnetar 
models with dipolar magnetic fields. As sketched in |Gabler et al.| 
( |2011| ) we extend this model by including an extended solid crust, 
which is able to support shear stresses. Here we complete the cor- 
responding discussion with a larger set of numerical simulations 
leading to quantitative estimates for the magnetic field strength. 

Our results reveal the existence of different regimes depend- 
ing on the magnetic field strength. For weak magnetic fields B < 
5 x 10 13 G we are able to recover purely crustal shear oscilla- 
tions. At intermediate fields 5 x 10 13 G< B < 10 15 G crustal 
modes are damped efficiently and the predominantly Alfven QPOs 
are limited to the fluid core of the neutron star. Finally, for strong 
fields B > 10 15 G, the magneto-elastic oscillations are no longer 
confined to the core by the interaction with the solid crust. They 
reach the surface with significant amplitudes and may be responsi- 
ble for the observed QPOs. In the analysis presented below we will 
make substantial use of an extension of the semi-analytic model 
presented in Cerda-Duran et al. (2009), which is based on a short 
wave-length approximation. We will further study the influence of 
the EoS on the properties of the equilibrium models, such as the 
change in frequency of the Alfven continuum or the dependence of 
the transition between the different regimes. 

The paper is organized as follows: in Section [2] we derive the 
equations for general-relativistic, magneto-hydrodynamics of elas- 
tic objects in the 3+1 split of general relativity. Those equations 
are then simplified to the case of torsional oscillations in the small- 
amplitude limit. Furthermore we discuss the boundary conditions 
we use, and the extension of the semi-analytic model of Cerda- 
Dura n et aT] ([2009). Section[3]is concerned with the equilibrium 
models we use as initial data for our simulations, and in Section 4 
we present the numerical methods employed. Section|5]contains the 
results of our simulation and we discuss our findings in Section[6] 
Finally, the appendices are devoted to derive an alternative method 
to compare our results with and to derive the eigenvalue problem 
of the shear modes in our model. 

In equations, we set c = G = 1, with c and G being the speed 
of light and Newton's gravitational constant, respectively. Latin in- 
dices run over (r, ip) and Greek indices over (£, r, ip). Partial 
and covariant derivatives are abbreviated by comma and semicolon, 
respectively. 



2 THEORETICAL FRAMEWORK 

In this section we present the equations governing the evolu- 
tion of torsional shear oscillations of magnetized neutron stars. 
First the equations of general-relativistic, magneto-hydrodynamics 
(GRMHD) are revised and the effects of an elastic crust up to linear 
order in the perturbations are included. Next, we apply a number of 
simplifications for small- amplitude, axisymmetric torsional oscil- 
lations. It turns out that these simplifications are equivalent to the 
so-called anelastic approximation, see Cerda-Duran et al. (2009). 
This provides a natural framework that can be used to generalize 
the formalism in use. The next subsection is concerned with the 
presentation of the boundary conditions applied in this work, and 
finally, we sketch the derivation of a semi-analytic model to calcu- 
late the Alfven continuum. 

Our numerical approach to study magneto-elastic oscillations 
of neutron stars is based on an approximate Riemann solver. Ad- 
ditionally, we have obtained a second method for the evolution of 
the crust coupled to the Riemann solver in the core. In this alter- 



native method we linearize the governing equations and evolve the 
resulting wave equation numerically. The results are used as a test 
to which we can compare our evolution with the Riemann solver 
method. We give the technical details of this approach in Appendix 

El 



2.1 General relativistic magneto-hydrodynamics and 
elasticity 

We adopt the framework of 3+1 split of general relativity. The gen- 
eral metric can thus be described by the following line element 



ds 



-adt 2 + ^ j {dx i + /3 i dt)(dx j + /3 j dt) , 



(1) 



where a and ft % are called lapse function and shift vector, respec- 
tively. For spherically symmetric space-times in isotropic coordi- 
nates the 3 -metric 7^ is conformally flat and can be related to the 
spatial flat 3 -metric 7^ by the conformal factor 

(2) 



Hi 



The description of matter in general relativity is based on the stress- 
energy tensor . For a magnetized perfect fluid with elastic prop- 
erties, can be written as a sum of different contributions: 



(3) 



(4) 



1 ~ 1 fluid ~r~ 1 mag ~T 1 e \ as • 

The different terms can be expressed as follows 

TZid = phuV + Pg»» , 

where p is the rest-mass density, h = 1 + e + P/p the specific en- 
thalpy, e the specific internal energy, P the isotropic fluid pressure, 
the 4-velocity of the fluid, and the metric tensor. Moreover, 
in the ideal MHD approximation, which we adopt here, the mag- 
netic part is given by 



mag 



;2 Lb V | 1 7 2 LbV 

b vru + -o g 



(5) 



with b^ 1 being the magnetic 4-vector and b 2 

The theory of elasticity in general relativity is based on 
the fundamental work of [Carter & Quintana (1972) and was re- 
cently extended in a series of papers by Karlovini & Samuelsson 
(2003 ] [20041 [2007]) and[Karlovini et al.| ( |2004| ). Our approach fol- 
lows Carter & Samuelsson H2006) . At the linear level in the pertur- 
bations the elastic effects of a crust enter in the form of the follow- 
ing contribution to the stress-energy tensor 



rpLtV 

elas 



(6) 



Here we have introduced the shear modulus us and the shear tensor 
E Ml/ . The Lie derivative of the latter along the 4-velocity £„E Ml/ is 
related to the strain tensor 



which is defined as 
where 

h"" := tf v + v?u v 



(7) 



(8) 



(9) 



is a projection perpendicular to the 4-velocity. 

The conservation of energy and momentum V 1 V T^ V — 
and the baryon number conservation V u J v — 0, together with 
Maxwell's equations V„ = and V U F^ = 4ttJ^ (when 
expressed in the ideal MHD approximation) lead to the following 
form of a flux-conservative hyperbolic system of equations 
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+ 



-g \ dt dx i 



(10) 



The quantities just introduced are the Faraday electromagnetic ten- 
sor field ', and its dual !F Ml/ , the rest-mass current := pu^ , 
and the electric 4-current J r/X , while yf^g — a- N /7 and 7 = 
det(7ij). The state vector U, the flux vector F\ and the sources 
S are given by 



U 



F l 



[D,Sj,T,B k ] : 



SjV l + 8) (P+ 



w 



v l B k -v k B l 



2 dx? 



aT» v 



<91na 



ab t B i 



,0,0,0 



(11) 



(12) 



(13) 



where v l = v l — /3 l /a and W = au l is the Lorentz factor. Addi- 
tionally, we introduce the 3-velocity of the fluid v l , the generalized 
rest mass density D, the momentum density Si, and the energy 
density r 



D — pW , 

Si = (ph + b 2 )W 2 v l 

r = (ph + b 2 )W 2 - 



- abib 1 , 



D. 



(14) 
(15) 

(16) 



The relation between B l , the magnetic field measured by an Eule- 
rian observer, and the magnetic 4-vector 6 M is given by 



WB l vi B t + W 2 v j B j v i 



W 



(17) 



The fluxes defined in Eq. |72| depend on the shear tensor 1T V 
which itself is not a function of the conserved variables U. To com- 
plete the system |T0| it is thus necessary to describe how the com- 
ponents of the shear tensor evolve with time. For simplicity, we will 
postpone the discussion of these terms to the next section, where we 
first apply some simplifications. 



2.2 Torsional oscillations in the small-amplitude limit and 
the anelastic approximation 

This work is concerned with torsional oscillations of neutron stars. 
Therefore we are allowed to perform a number of simplifications to 
the full set of ideal GRMHD equations presented above. 

(i) The temperature of an old neutron star is well below the 
Fermi temperature, i.e. a barotropic EoS describes the matter with 
sufficient accuracy. Hence, the equation for the energy density r in 
system (To| becomes redundant. 

(ii) For purely poloidal background fields and in axisymmetry, 
the torsional oscillations decouple at the linear level from the polar 
ones. For small- amplitude oscillations it is thus justified to evolve 
B^ and S^ separately, and to keep B r , B e , S r , and Sg constant in 
time. 

(iii) When setting the initial velocities in the r and 6 directions, 
and therefore the polar perturbations, to zero, they remain zero dur- 
ing the evolution (see (ii)). The continuity equation reduces to 



dt 



dt 



= 0, 



(18) 



because we assume axisymmetry, i.e. the density D stays constant 
throughout the evolution. 



(iv) Because we are dealing with torsional oscillations, which 
do not involve large density changes, see (iii), and which only cou- 
ple very weakly to the spacetime evolution, the Cowling approxi- 
mation should hold, i.e. the metric can be regarded as fixed. Fur- 
thermore, the matter distribution and thus the spacetime can be re- 
garded as spherically symmetric, since the deviations from spheric- 
ity due to the magnetic field are small for the magnetic fields con- 
sidered here ( Bocquet et al. 1995] see). In the case of isotropic co- 
ordinates and for non-rotating stars this implies that the shift vector 
vanishes, ft 1 — 0, and consequently v z — v % . 

In general, the eigenvalues of the flux-vector Jacobians of the 
system |T0] ) depend on the metric tensor, the conserved variables, 
and the speed of sound. The latter dependence strongly limits the 
numerical time step, slowing down any time-explicit numerical 
simulation. To overcome this limitation in a more general context, 
Bonazzola et al. (2007 ) suggested to remove all pressure dependen- 
cies from the fluxes F l , which is known as the anelastic approxima- 
tion. This changes the properties of the eigenvalues of the system 
making them independent of the speed of sound. In our case, as a 
consequence of (iii), the eigenvalues of the flux- vector Jacobians 
of the system |T0| do no longer depend on the speed of sound (see 
below in this section) and the numerical time- step is thus set by the 
speed of Alfven waves, only. 

The anelastic approximation provides a natural generalization 
of the method applied in the present work. It allows one to describe 
situations also involving perturbations in density and pressure, and 
to study more complicated magnetic field configurations including 
the coupling between axial and polar oscillations, as a first step 
towards full 3-dimensional simulations. This is, however, beyond 
the scope of the present work. 

To summarize, we only evolve the equations for S^ and B v 
and hold all other conserved variables of Eq. (IT) constant. To- 
gether with Eq. |T0| the following definitions describe the evolu- 
tion of torsional oscillations of magnetized neutron stars including 
an elastic crust: 



U 
F r 



[S<p,B*], 
b^B r 
W 



><pB° 
W 



2/isS r 



■ 2/xsS" 



-v*B r 



-v*B b 



[0,0]. 



(19) 
(20) 

(21) 
(22) 



To derive the shear tensor for torsional oscillations in 
the small- amplitude limit, we follow Schumaker & Thorne (1983). 
Since we assume the Cowling approximation to hold, a diagonal 
metric g^ v — diag^^), and restrict our considerations to linear 
order in the displacements, the components of Y>^ v can be written 
as 



_L_ 

2^ 



+« J *(«'«)J -£(<*») 



(23) 



while the other components vanish, E Mt = Y? v — 0. Here, the dis- 
placement ^ has been introduced. The corresponding time deriva- 
tive is related to the 3-velocity of the fluid as follows 



(24) 



For purely torsional oscillations the more general expression 
in Eq. ([23]) simplifies to 
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g rr ^ r 

g 66 e,e 

9 rr ^ r 9 ee e e 



(25) 



In its current form the system of equations |T0| and fT9| - 
(22} is incomplete, because the displacement appears in the fluxes, 
Eq. (20} and $2l\ . Hence, it is necessary to add (at least) one equa- 
tion for the evolution of During the calculations it turned out 
that the direct integration of £^ via Eq. (24} is numerically disad- 
vantageous because of several reasons: 

(i) The use of requires an evaluation of its spatial deriva- 
tives when calculating the fluxes in fT9} - (22} . This leads to nu- 
merical inaccuracies in regions where the displacement changes 
strongly. 

(ii) Alternatively one could move the additional terms includ- 
ing the derivatives of ^ to the source terms of the system of equa- 
tions. However, it is numerically more desirable to have a system of 
conservation laws without any source terms. Furthermore, the spa- 
tial derivatives of the shear modulus /is would appear in the source 
terms, which might give rise to spurious oscillations, because /is is 
given in tabulated form. 

(iii) The interface conditions at the crust-core interface are 
given in terms of and £^ r . It turns out that it is numerically 
difficult to impose them directly on ^ ensuring a sufficiently ac- 
curate corresponding condition for its radial derivative (see end of 
Section[4}. 

Therefore, instead of evolving ^ (p , the system of Eqs. \\0\ with 
Eqs. (T9] > - ( |22| is supplemented with the following equations for 
the spatial derivatives of ^ 



(ex* 



0, 
0, 



(26) 
(27) 



which are based on Eq. (24} . 

To obtain the eigenvalues we write the complete system in the 
form of a conservation law like Eq. (TO}. This introduces the new 
conserved variables U, fluxes F k , and sources S: 



U 



F = 



(28) 



(29) 



(30) 



where k = {r, 0}. The solution of the eigenvalue problem of the as- 
sociated flux- vector Jacobian leads to the following non-zero eigen- 
values 




M/2 



A 



A 

: phW A (l + vyv*) + B r B r + B e B e 



(31) 



(32) 



The set of equations (10} with (28}, (29}, and (30} is useful for 
computing the eigenvalues, but contains non-zero sources in the 
conservation law. For the numerical time-evolutions, we thus im- 
plement the original set Eq. (To} with Eqs. (19} - (22} and Eqs. (26} 
and (27} (see Section[4]below for details). 



2.3 Boundary conditions at the surface and treatment of the 
crust-core interface 

The system of equations for magneto-elastic torsional oscillations 
contains two degrees of freedom related to the two non- vanishing 
eigenvalues (31}. Therefore, we have to impose boundary condi- 
tions at the surface of the star that mimic the incoming waves from 
the magnetosphere, which are not included in our simulations. In 
the core, where /is = 0, there are still two degrees of freedom. 
Since we are simulating both regions (crust and core) having the 
same number of degrees of freedom there is no need for boundary 
conditions at the crust-core interface (although there is a need for a 
special treatment for numerical reasons, see below). 

In addition, in the case of ideal MHD without charges, the 
electric field is continuous everywhere, and hence the velocity v^, 
too. This implies continuity of the displacement, and of its time 
derivative, £^ t . At the surface of the star and at the crust-core inter- 
face the tangential derivative, is continuous, while no restric- 
tions apply to the continuity of £f r . 

Boundary conditions at the surface: We assume that there are 
no current sheets at the surface of the star, i.e. the tangential mag- 
netic field components have to be continuous 



(33) 



^crust ^atmosphere 

at the surface. 

The conservation of momentum gives the continuity of the 
traction 



T(h, (p) — T(f , (p) — T rif 



(34) 



i.e. the tangential stresses inside and outside the star have to balance 
each other. Here the tilde indicates normalized vectors, n is the 
normal to the surface of the star and thus n = f . 

The continuous traction condition can be simplified in the case 
of continuous W and leads to a condition for £^ r : 



u ^crust 



+ 



Ms, 



rrinp 

- 1 crust 


atmosphere 


(35) 


crust, r 


— b r b v 

atmosphere 


(36) 


crust, r 


= 0. 


(37) 



Eqs. (33} and (37} are the set of boundary conditions that we apply 
at the surface of the star. We need the additional condition (37}, be- 
cause we are evolving more variables than the degrees of freedom 
of the system. 

The assumptions made here are motivated by the picture that 
the magnetospheric field close to the surface will move with its 
footpoints in the crust, i.e. the exterior solution relaxes to a force- 
free field on a much shorter time scale than the interior evolves. 
This implies that currents can be maintained in the magnetosphere, 
which is necessary to support more general equilibrium configu- 
rations than considered here and hence to create a twisted magne- 
tospheric field. A more detailed discussion of the coupling to the 
magnetosphere would exceed the scope of this paper and will be 
purpose of further investigations. 

Our boundary conditions are similar to those used in previous 
work without the presence of a crust ([Sotani et a l. 008b, Cerda- 
|Duran et a l. 2009 ; Colaiuda et al. 2009 ) and in simulations with a 
crust ( |Gabler et aEl|2011| |Colaiuda & Kokkotas||2011| ). However, 
our boundary conditions differ from those of Lander et al. ( 2010 ) 
and |Lander & Jones (2011} who choose a particular set of vari- 
ables vanishing at the surface of the neutron star, which defines the 
corresponding boundary conditons. Other approaches ( Braithwaite 
|& Nordlund|2006||Lasky et al.|2011||Ciblfi et al.|2011| involve the 
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evolution of some parts of the neutron star's atmosphere. For purely 
toroidal oscillations it is possible, however, to impose appropriate 
boundary conditions and to avoid that evolution. 

The treatment of the crust-core interface requires particular 
attention. At this interface no boundary conditions are required and 
by knowing the variables at one instant of time on both sides (crust 
and core) one should be able to evolve the system. However, the 
stability of the employed schemes turned out to depend sensitively 
on the particular treatment of the reconstruction of the variables, 
which are allowed to be discontinous or to have discontinous spatial 
derivatives (see below). In the following we will describe how to 
proceed to achieve stable evolutions. 

As anywhere else the traction is supposed to be continous 



d 2 Y( X ,C,t) 



rprtp 
core 



rpr^p 
crust : 



u u crust ^4 S ,r ■ 



(38) 
(39) 



This can be transformed by virtue of the linearized induction equa- 
tion 

b* = b r Z* r + b 9 ^ e , (40) 
the continuity of the displacement, and thus Ccore e ~ £ crust 6" mt0 



' ^4 > crust, r i 



core.r 



1 + 



Ms 



Ql(pr\2 I S crust.v 



(41) 



(42) 



In Section [4] we will show how to ensure the continuity of the trac- 
tion and maintain the correct relation between the radial derivatives 
of the displacement in the core and in the crust. 

Obviously the discontinuous radial derivative allows W to be 
discontinuous (Eq.[40j. Hence, in general, there are current sheets 
present at the crust-core interface, which are unavoidable, and a 
consequence of the coupled evolution and the assumption of ideal 
MHD. 



&Y(x,C,t) = 

o>t 2 ttot(x) 2 dC 2 

Next we assume standing waves of the form 
y (X, C t) = a(x) sin« + 4>c) cos(2tv ft + 00 : 



(44) 



(45) 



where k is the wavenumber, a(x) the amplitude, <j> t the temporal 
phase, a spatial phase, and / is the oscillation frequency. The 
dispersion relation then is simply given by 



/ 



(46) 



At this point the frequencies of the oscillations are completely de- 
termined by the magnetic field topology and the boundary condi- 
tions. 

Cerda-Duran et al. (2009 ) did not consider an extended crust, 
and hence the boundary condition was set at the surface and corre- 
sponded to the continuous traction condition. This resulted in a van- 
ishing radial derivative of the displacement and a maximum ampli- 
tude of the perturbation at the surface. However, when an extended 
crust is present there exist two different regimes. For low magnetic 
field strengths (B < 10 15 G) the standing waves show a node at the 
crust-core interface ( Gabler et al. 2011), which may be interpreted 
as a reflection of the standing wave at the crust-core interface. As 
we show below, this change of the boundary condition for the semi- 
analytic model is necessary to calculate the correct frequencies and 
to find the symmetry of the numerically obtained QPOs. 

For stronger magnetic fields (B > 10 15 G), the oscillations 
reach the surface, and we can apply the boundary condition of 
Cerda-Duran et al. ( 2009). In order to take the crust into account 
the velocity of the perturbation can be approximated by the eigen- 
values ( [37] ), assuming that the perturbation Y is still traveling along 
the magnetic field lines. In an intermediate regime at around 10 15 G 
we do not expect the semi-analytic model to be valid, because in 
this case the shear and magnetic contributions to the evolution in 
the crust are of similar order. 



2.4 Semi-analytic model 

A very useful tool to obtain the frequencies of purely magnetic os- 
cillations inside a neutron star was presented by Cerda-Duran et al.| 
(2009 ). These authors showed that in the linear regime and in the 
limit of short wavelengths it is possible to calculate the Alfven con- 
tinuum with a semi-analytic model. Here, we will only sketch the 
method, and for more information we refer to Cerda-Duran et al.| 
(2009 ). In the aforementioned limit an Alfven wave will travel 
along magnetic field lines corresponding to 

£=v.(x), (43) 

where v a is the Alfven velocity. Any displacement Y traveling 
along the magnetic field lines can be expressed as a function of 
magnetic -field-line -adapted coordinates (%, Q and time, where x 
labels the magnetic field line by the radius at which it crosses the 
equatorial plane, and f = £(r, 0; x) Atot (x) — 1/2 is a dimension- 
less parameter along each field line. Here t tot (x) is twice the total 
travel time of an Alfven wave traveling along a magnetic field line 
starting from the equatorial plane and ending at the surface or at 
another point in the equatorial plane. Note that ( used in this work 
corresponds to £ in |Cerda-Duran et aT| (p009), because here £ de- 
notes the displacement related to the 4-velocity of the fluid. For a 
traveling wave Y satisfies trivially the wave equation 



3 EQUILIBRIUM MODELS 

The initial models are self-consistent general relativistic equilib- 
rium models of magnetized non-rotating neutron stars with a purely 
poloidal magnetic field ( Bocquet et al. 1995). We use the numerical 
code "magstar" of the LORENE librar^to compute these models, 
which include the effects of the magnetic field on the matter and 
the spacetime. 

In our models the magnetic field is generated by a current of 
the form — phC, where C is a constant which determines 
the strength of the magnetic field. We note that different choices 
of the current can lead to different internal magnetic fields, still 
resembling an exterior dipole field (see e.g. van Hove n & Levinj 
2011). Hereafter, we will label the different models by the surface 
value of their magnetic field strength at the pole. 

As we assume a spherically symmetric spacetime and matter 
background in our simulations we angle-average the density of the 
background model to obtain a spherically symmetric model from 
the LORENE data. In the most extreme cases, i.e. for a very strong 
magnetic field this simplification changes the structure of the neu- 
tron star by at most about one per cent. For example the density at 
different angles but constant radius varies within less than one per 

1 http://www.lorene.obspm.fr 
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EoS 


mass 


circunferencial 


inner radius 


relative size 






radius r s 


of crust 


of crust 




[M ] 


[km] 


[km] 


Ar crust j-o^j 


APR+DH 


1.4 


12.10 


11.22 


7.2 




1.6 


12.07 


11.38 


5.7 




1.8 


12.00 


11.43 


4.8 




2.0 


11.90 


11.44 


3.9 




2.2 


11.63 


11.31 


2.8 


APR+NV 


1.4 


11.94 


10.85 


9.1 




1.6 


11.93 


11.07 


7.2 




1.8 


11.92 


11.19 


6.1 




2.0 


11.81 


11.23 


4.9 




2.2 


11.56 


11.12 


3.8 


L+DH 


1.4 


14.74 


13.33 


9.6 




1.6 


14.85 


13.72 


7.6 




1.8 


14.93 


13.94 


6.6 




2.0 


14.99 


14.13 


5.7 




2.2 


14.94 


14.22 


4.8 


L+NV 


1.4 


13.29 


11.88 


10.6 




1.6 


13.58 


12.35 


9.1 




1.8 


13.86 


12.76 


7.8 




2.0 


14.02 


13.08 


6.7 




2.2 


14.12 


13.30 


5.8 



Table 1. EoS, masses, radii of the star, radii of the crust-core interface, and 
sizes of the crust of the models studied in this paper (without magnetic 
field). 



cent at B — 5xl0 15 G. Therefore, the influence on the oscillations 
is in general less than that of other approximations we are applying. 

For the EoS we can choose between different realistic 
barotropic models including the description of a crust. We use four 
combinations of two EoS in the core matched to two distinct EoS 
for the crust. For the core we chose the APR EoS ( |Akmal et al.| 
[1998) and the stiffer EoS L jPandharipande & Smith|T975|), while 
for the low density region of the crust we select EoS NV ([N egele 
|& VautherinfL973] ) and EoS DH ( |Douchin & Haensel|200l])7Th e 
recent discovery of a 2 M© neutron star by Demorest et al. (2010), 
excludes EoS which cannot reproduce such large masses. The prop- 
erties of the equilibrium models used in this work are summarized 
in Table|T] These models are a subset of the models used in Sotani 
|et al.|p007] >. Here, r s and r cc are the radii of the surface of the star 
and of the crust-core interface, respectively, and Ar cr ust = r s — r cc 
is the size of the crust. We note that table 1 of Sotani et al. (2007) 
shows Ar C rust /rcc, instead of Ar cr ust/r s and thus the percentage 
for the relative size of the crust is different in their case. However, 
we checked that the value of Ar cr ust is the same in both cases. The 
frequencies of the crustal modes with n > depend sensitively on 
the size of the crust (Samuelsson & Andersson 2007 ). Therefore, a 
proper definition of the size of the crust is important. 

The shear oscillations of the crust are mainly determined by 
the shear modulus /is which we obtain from the zero-temperature 
limit of Strohmayer et al. (1991) given by 



p s = 0.1194 



n x (Ze) 2 



(47) 



where m is the ion density, (Ze) the ion charge and a = 
[3/(47rni)] 1 ^ 3 the average ion spacing. This equation is derived 
for a perfect bcc lattice, and the shear modulus, which has different 
magnitude along different crystal axes, is averaged in order to ob- 
tain an isotropic effective shear modulus /is (see Strohmayer et al.| 
|1991| >. For the NV EoS of the crust we use a simple fitting formula 
derived by |Duncan|(T998] > 



/is = 1-267 x 10 dO ergcm" d pIi 5 , (48) 

where p±4 = />/(10 14 g cm -3 ). 

To calculate the shear modulus for the DH EoS, one has to 
evaluate m in Eq. (47J in terms of the nucleon number A, the proton 
number Z, and the neutron fraction X n (Piro 2005): m = pivra 
and Api ~ p(l — X n ). The composition at a given density is given 
in Douchin & Haensel (2001). The shear modulus can be estimated 
to be 



/is = 1.2 x 10 erg cm 



X n 



0.25 



-3 4/3 
Pl4 



4/3 



/302\ 



38 ) { A ) 



(49) 



|Sotani et al. (2007) introduced the following fit to this equation 

Ms 



1 r>30 

10 erg cm 



3 (0.02123 + 0.37631pi4 + 3.13044/>? 4 



4.718141p? 4 + 2.46792/4) 



(50) 



This function provides a good approximation for densities larger 
than p = 5 x 10 n gcm -3 , but below we will rely on the more 
general expression in Eq. \49\ . The main motivation to use this fit 
is to allow for a direct comparison of the results obtained in this 
work to the results presented in |Sotani et al.| ( [2007] ) and |Col aiuda 
|& Kokkotas| ( [20Tl] ). 

The crust-core boundary for the NV and DH EoS is defined at 
p CC)NV = 2.4 x 10 14 gcm- 3 , and/>cc,DH = 1.28 x 10 14 gcm~ 3 , 
respectively. 

We employed two EoS which give neutron star models with 
large shear moduli compared to other available EoS ( Stei ner &| 
Watts 2009 ). Therefore, the results obtained in this work can be 
regarded as an upper limit for the influence of the crust. 



4 NUMERICAL METHODS 

This section is devoted to the numerical methods employed to ana- 
lyze the torsional shear oscillations. The current work is an exten- 
sion of the study of Cerda-Du ran et al.| ( |2009| ), whose non-linear 
GRMHD code we use as a basis for our simulations. The code has 
been developed in order to investigate various astrophysical scenar- 
ios where both magnetic fields and strong gravitational fields play 
an important role in the evolution of the system. ( [Dim melmeier 
|et al.|2002a|bl[2005l|Cerda-Duran et al.|2008) . The code uses high- 
resolution shock-capturing schemes to solve the GRMHD equa- 
tions for a dynamical spacetime, under the approximation of the 
conformally flat condition (CFC) for the Einstein's equations jlsen-| 
|berg|2 008 ; Wil son et al.|1 996). For a spherically symmetric space- 
time the CFC metric is an exact solution of Einstein's equations, 
and reduces to the solution in isotropic coordinates. Therefore, this 
numerical code is well suited to describe the spacetime used in the 
simulations of the present work. The equations are cast in a first- 
order, flux-conservative hyperbolic form, supplemented by the flux 
constraint transport method to ensure the solenoidal condition of 
the magnetic field. 

The basic version of the code including the solution of the 
ideal GRMHD equations was thoroughly tested in Cerda-Duran 
|et al.| ( |2008| ), who demonstrate the robustness of the code for a 
number of stringent tests, such as relativistic shocks, highly mag- 
netized fluids, equilibrium configurations of magnetized neutron 
stars, and the magneto-rotational core collapse of a realistic pro- 
genitor. One important feature is the ability of the code to handle 
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different classes of EoS which range from simple analytical expres- 
sions to microphysically derived tables. We want to emphasize that 
although the current project is concerned with small-amplitude per- 
turbations in order to apply simplifications appropriate to a linear 
regime, the code can in principle handle large amplitudes and in 
general is nonlinear. 

While the numerical method employed in the original version 
of the code to integrate the equations is unchanged, we modify 
the equations to be solved according to Eq. p0| and {T9}-(|22} in 
the crust. As mentioned above we also have to solve for the spa- 
tial derivatives of ^ <p . For these equations, which do not represent 
conservation laws, we calculate the fluxes — av^ corresponding to 
Eqs. ( [26} and (27} at the cell interfaces with the corresponding ap- 
proximate Riemann solver. The derivatives of the fluxes are ap- 
proximated by dividing the difference of the flux at the two cell 
interfaces by the grid extension. Although this approach should be 
able to cope with discontinuities, as present at the crust-core inter- 
face for the shear modulus, it turns out that special care has to be 
taken to achieve a converging method. The coupling between the 
crust and the core will be described in detail below. We note that 
when using this approach and setting the shear modulus in the crust 
to zero the results of Cerda-Duran et al. ( 2009) are recovered. 

Setting up the interface conditions turned out to be a very deli- 
cate issue. The shear modulus is discontinuous at the crust-core in- 
terface, and Riemann solvers are able to cope with discontinuities 
at cell interfaces. Therefore, we define the crust-core interface to 
be located at a cell interface. In this case it is crucial to ensure that 
the reconstruction procedure gives a value for £f r ~ E 7 ^ which 
is consistent with the continous traction conditions (Sec. |2.3| ). Any 
standard reconstruction method not taking this condition explicitly 
into account failed and the simulations produced spikes in the radial 
profiles that spoiled the evolution or even made the whole evolution 
unstable. The main cause for this behavior is the discontinuity of 
the shear modulus at the crust-core interface. For intermediate and 
weak magnetic fields, the very large shear modulus on one side and 
the vanishing shear modulus on the other causes the different terms 
in the momentum equation for the radial flux Eq. ( [20} to be much 
larger on the side of the crust than on the side of the core, i.e. 



b v B r 



W 



Ei r 



V I crust 



>> 



b^B r 



W 



(51) 



The evaluation of the flux at the crust-core interface is numerically 
problematic due to non-cancellations of the two terms on the side 
of the crust. However, when taking the continous traction condition 
appropriately into account, this problem does not arise and the flux 
at the core-crust interface is well behaved, allowing one to perform 
simulations also for intermediate and weak magnetic fields. 

We are using the following numerical treatment based on the 
continous traction condition ^ ore = £ crust = ^ an ^ £core,r — 
r]^ c rus t, r with ?7 = 1 + /is / (& 4 (b r ) 2 ). The derivatives at the crust- 
core interface can be approximated by 

^(^cc) — £core(^C' 



S cr 



0.5Ar 
+ 0.5Ar) 



0.5Ar) 

i 

-r(fcc) 



(52) 



(53) 



£core,r (^cc) 
f crust.r^cc) ().5Ar 

where Ar the grid spacing in radial direction. Both equations lead 
to the following expression for 

_ £Sore(rcc - 0.5Ar) + 7^ rust (r cc + 0.5Ar) 



l+ri 



(54) 



Knowing ^ at the crust-core interface one can calculate the ra- 



dial derivatives £ 



crust, r 



and £core,r> and finally the fluxes. For the 



calculations presented in this paper we used a second-order approx- 
imation of the derivatives instead of Eqs. ([52} and ([53}. 



5 RESULTS 

First, we show that our code reproduces the purely shear oscilla- 
tions of the crust, and we demonstrate how these crustal modes 
are damped when the magnetic field is switched on. The next two 
subsections are concerned with the behavior of the magneto-elastic 
oscillations at intermediate fields (5 x 10 13 G< B < 10 15 G), and 
at higher magnetic field strengths, respectively. In the remaining 
subsections we address the issue of different equilibrium models 
and investigate the magnetic field threshold at which the magneto- 
elastic QPOs break out of the core and reach the surface of the 
neutron star. For intermediate magnetic field strengths the QPOs 
are confined to the core due to the interaction with the crust. If not 
stated otherwise we refer to the results using an equilibrium model 
obtained with the APR+DH EoS with a mass of M = 1 .4 M (see 
Table[T}. 

The 2-dimensional simulations are computationally very de- 
manding. On a single Intel Xeon processor with 3GHz a typical 
simulation necessary to perform a Fourier analysis of the Alfven 
spectrum takes of the order of a few weeks for one equilibrium 
model. We are therefore restricted in the maximum time we can 
evolve the models. In the most demanding situation the minimum 
evolution time corresponds to 10 times the dynamical time scale 
which is set by the Alfven crossing time. For a magnetic field of 
10 14 G the Alfven crossing time for the fundamental Alfven oscil- 
lation is about 2 s. This value scales inversely proportional to the 
magnetic field. 



5.1 Recovering purely shear oscillations of the crust 

The purely shear oscillations for various realistic EoS in gen- 
eral relativity have been calculated for the linearized problem 
by |Messios et al.| ( [2001} and |SotanT"e t al. (2007). To recover 
their results, we performed a series of simulations for a selec- 
tion of models with zero magnetic field strength. As initial ve- 
locity perturbation we use a simple radial law in the form v w 
sin (tt/2 * (r — r cc )/(r s — r cc )) multiplied by a sum of the first 
ten vector spherical harmonics for the angular dependence. With 
this kind of perturbation we ensure to excite single modes with dif- 
ferent values for the radial and angular mode numbers n and /. We 
evolve the system for 1 s in the case of n — modes and 50 ms 
for the modes with n ^ 1. The resolution of the simulations was 
120(r) x 60(0) points in the domain [0, r s ] x [0,7r/2], the grid 
is equidistant in both directions and we used equatorial symmetry. 
The chosen grid corresponds to about 20 radial zones inside the 
crust. 

Table[2]gives the oscillation frequencies of our dynamical sim- 
ulations extracted from the Fourier analysis at points inside the 
crust. The frequencies of the modes agree up to a few percent with 
those of the linear approximation (given in round parenthesis). For 
modes with n ^ 1, the frequency resolution of the Fourier trans- 
form does not allow us to resolve modes with the same n but dif- 
ferent /, which are only separated from each other by a few Hz. 
Therefore, the measured frequency is a mixture of different / con- 
tributions, and its value is expected to be slightly larger than that of 
Sotani et al. ( 2007) for / = 2, which are reported in the Table[2]as 
well. 
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Model mode frequency in Hz 



n=0(±lHz) n=l n =2 





1 = 2 


1 = 3 


1 = 4 


1 = 5 


1 = 6 


(±20Hz) 




APR+DH 1.4 
APR+DH 1.6 
APR+DH 2.0 
L+DH 1.6 
L+DH 2.0 


25.4 (24.6) [25.1] 
24.3 (23.4) [24.0] 
21.9 (21.3) [21.7] 
21.0 (20.6) [20.9] 

19.5 (18.9) [19.2] 


40.0 (38.9) 

38.5 (37.0) 

34.6 (33.6) 

33.1 (32.5) 

30.7 (29.9) 


53.6 (52.2) 

51.2 (49.6) 

46.3 (45.1) 

44.8 (43.7) 

40.9 (40.2) 


67.3 (65.1) 
64.3 (61.9) 
58.0 (56.3) 
55.6 (54.5) 
51.2 (50.1) 


80.0 (77.8) 
76.5 (73.9) 
69.2 (67.3) 
66.8 (65.1) 
61.4 (59.9) 


741 (761) [734] 
829 (860) [825] 
1052 (1083) [1045] 
565 (586) [567] 
682 (713) [677] 


1190(1270) 
1340(1430) 
1842(1810) 
917 (980) 
1100(1190) 


APR+NV 1.6 
APR+NV 2.0 
L+NV 1.6 
L+NV 2.0 


23.9 (23.8) [23.6] 
21.5 (21.4) [21.2] 
22.0 (21.8) [21.8] 
19.5 (19.7) [19.6] 


37.5 (37.6) 
33.7 (33.9) 
34.7 (34.5) 
31.2 (31.1) 


50.7 (50.5) 
45.4 (45.5) 

46.8 (46.3) 

41.9 (41.7) 


62.9 (63.0) 

56.5 (56.7) 

58.6 (57.7) 
52.2 (52.1) 


75.5 (75.3) 
67.8 (67.8) 
69.7 (69.0) 
62.5 (62.2) 


692 (689) [684] 
838 (858) [827] 
526 (525) [522] 
615 (615) [608] 


1230(1220) 
1501 (1520) 
936 (930) 
1092 (1090) 



Table 2. Frequencies of some torsional shear modes of the crust for different EoS. Numbers following the abbreviation of the EoS give the mass of the stellar 
model in solar mass units. The frequencies in round parenthesis are from Sotani et al. ( 2007 1, and the squared brackets give the result of the eigenmode analysis 
(see Appendix[B}. For the n = 1 and n = 2 modes we compare to the Z = 2 case only. The frequencies of the eigenvalue calculation for n = and / > 2 
can be obtained by multiplying the corresponding frequencies for I = 2 with y/ (I — + 2) / 2 . The error ranges shown in the table header correspond to 
the frequency spacing in the Fourier analysis. 



%jr 0.6 



q- ^0.4 
0.2 




- 90x30 

- 120x40 (c) 
150x50 (m) 

- 180x60 (f) 




400 600 
time [ms] 



1000 



Figure 2. Numerical damping of crustal oscillations. The upper panel shows 
the evolution of the maximal amplitude of £^ normalized to its initial value, 
at a point in the crust near the pole for different grid resolutions. The numer- 
ical damping decreases with increasing resolution. The lower panel shows 
the order of convergence (circles) computed using the three highest res- 
olution simulations compared to the expected second order convergence 
(dashed line). 



Additionally, we have computed the frequencies from the as- 
sociated eigenvalue problem (Appendix[B]>. For the computation of 
these frequencies we used a radial grid of about 80 zones in the 
crust. The results for the n = 0, / = 2 and n = 1 modes are given 
in squared brackets in Table[2] They agree up to a similar accuracy 
with those that have been obtained in the literature. 

To investigate the convergence properties of our code when 
no magnetic field is present, we have calculated the evolution of 
the n — and / = 2 mode of our reference neutron star model 
for different grid resolutions: 180 x 60, 150 x 50, 120 x 40, and 
90 x 30. The angular grid is equidistant, while the radial grid is 
equidistant only in the crust, where 40 percent of the zones are lo- 
cated, and coarsens towards the center of the star. The finer mesh 
in the crust ensures higher accuracy without significant increase 
of the computational costs. In order to save further computational 



power we assume equatorial symmetry. The mode frequencies ex- 
tracted from the simulations at different resolutions agree within 
the frequency resolution of the Fourier transform. The upper panel 
of Fig. [2] shows the time evolution of the maximum amplitude of 
at the crust for different grid resolutions. One clearly sees that 
the numerical damping decreases with increasing resolution. For 
this simple test case it is possible to compute the order of conver- 
gence using the results for the three highest resolutions, when one 
assumes that the error in the interesting variable, /, scales as A p , 
where A is the size of the numerical cell and p the order of conver- 
gence. To compute p one searches for roots of 



medium 



fn 



/fin 



medium 



A£ 

fine 



(55) 



where the subscripts denote fine, medium or coarse grid resolution. 
The lower panel of Fig. [2] shows that after a short initial transient 
the order of convergence, p, rapidly converges to 2, which is the 
expected order of convergence of our numerical scheme. 



5.2 Absorption of crustal shear modes by the Alfven 
continuum 

This subsection is concerned with the absorption of purely shear 
modes of the crust into the Alfven continuum of the core. We will 
often refer to this process as the damping of crustal shear modes. 
This expression may mislead the reader to think of dissipation pro- 
cesses. However, damping in the current context refers to the trans- 
fer of energy from crustal modes into the continuum of the core, 
but not to dissipation of energy. Only when referring to numerical 
damping we mean the usual concept describing the loss of ordered 
kinetic energy by numerical dissipation. 



5.2.1 n = shear modes 



As we have shown in Ga bler et aT] |201 1 ) purely shear oscillations 
are absorbed very efficiently by the Alfven continuum of the core 
for magnetic field strengths B > 5 x 10 13 G. In this case the am- 
plitude of the perturbations of the crust is damped by transferring 
their energy to the Alfven continuum. To analyze how this damping 
scales with the magnetic field strength, we have performed a series 
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200 300 
time [ms] 



500 



Figure 3. The maximum of the normalized overlap integral for the n = 0, 
/ = 2 eigenmode of the crust as a function of time. Stronger magnetic field 
results in faster damping of the initially excited crust mode. Orthogonality 
of different eigenmodes is fulfilled numerically up to about 10 -5 . In the 
legend we introduced the abbreviation Bis = 10 13 G. 



0, 



^ 2 for 

13 



of simulations for different crustal modes with n - 
different magnetic fields (0, 10 13 ,2 x 10 13 , 5 x 10 13 , 8 x 10 13 ,10 14 
and 2 x 10 14 G). For these simulations we use a grid of 150 x 100 
zones covering a domain [0, r s ] x [0, tt], which is equivalent to the 
medium grid of the previous section but not assuming equatorial 
symmetry. Here, we use the solution of the eigenvalue problem for 
the unmagnetized crust (see Appendix[B]> as initial perturbation for 
the velocity. Symmetries are exploited whenever a perturbation is 
purely symmetric or antisymmetric with respect to the equatorial 
plane. 

In the following we will discuss our results by using so-called 
overlap integrals (derived in Appendix[B]>. These overlap integrals 
are the expansion coefficients of an arbitrary spatial function in the 
basis of the crustal oscillation eigenmodes, i.e. they give a measure 
of how strong the different eigenmodes are excited (see Gable r et al.| 
2009 for an application to radial oscillations of neutron stars). In 
Fig.[3]we show the time dependence of the maxima of the overlap 
integrals defined in Eq. JbTTJ corresponding to the n = 0, I — 2 



crustal mode for different magnetic field strengths. The stronger 
the field is the faster the damping of the shear mode proceeds. For 
high field strength, B > 2 x 10 14 G, it is not possible to obtain 
a characteristic damping time r, because the time scale is shorter 
than one oscillation period. Therefore, the latter can be used as an 
upper bound for r. For 2 x 10 14 Gwe show the evolution of the 
overlap integral only up to the time when global magneto-elastic 
oscillations start to dominate and interfere with the purely shear 
modes of the crust. 

Table[3] shows the damping timescales, obtained by analyzing 
the overlap integrals, for different n = 0, / ^ 2 modes and different 
magnetic field strengths. The values for zero magnetic field serve 
as a measure of the numerical damping of the code. As expected 
for numerical dissipation processes, modes with higher / suffer 
stronger from numerical damping than lower / (see also [Cerda- 
|Duran|2oTo] ). For weak magnetic fields, e.g. 10 13 G, the damping 
of high I modes is dominated by numerical dissipation, while for 
low / modes it is caused by the interaction with the Alfven contin- 
uum of the core. For magnetic fields stronger than 5 x 10 13 G, we 
are confident that the damping time of all studied modes is physical 
and not due to numerical dissipation. Above 2 x 10 14 G the oscilla- 
tions are damped on shorter time scales than the respective oscilla- 
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Table 3. Damping timescale r in ms for different n = 0, crustal shear 
modes for different magnetic field strengths. 
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Figure 4. The dependence of t/£a on the magnetic field. The stronger the 
magnetic field is the less is the spread of the numerical values for differ- 
ent / around this value. The shaded area indicates the maximum range of 
variation of t/£a around the mean value (~ 0.04). 



tion period, and hence it is impossible to obtain accurate damping 
times. 

Figure [4] shows t/£a, where £a is the Alfven crossing time 
of the star at the pole. The damping time of crustal modes due to 
the absorption by the Alfven continuum scales linearly with £a, 
i.e. r decreases with increasing magnetic field. The mean damp- 
ing time is about 0.04 £a- Deviations from this value (see Fig.[4} 
depend non-trivially on the magnetic field and the mode number I. 
The spread decreases with increasing magnetic field strength, being 
smallest for our B = 2 x 10 14 G simulation. Low /-modes (filled 
circles in Figure[4]) show a smaller spread around the mean value 
than the higher /-modes (crosses in the same figure). 

These deviations are expected, because the damping depends 
on a variety of parameters as for example the frequency of the 
crustal mode, the frequencies available in the Alfven continuum of 
the core, and the spatial structure of the crustal modes. Numerical 
effects may also affect the damping times. These are: Firstly, the 
grid resolution necessary to obtain a comparable accuracy for dif- 
ferent modes increases with increasing mode number /. Secondly, 
waves in the crust reaching the core-crust interface will propagate 
into the core as Alfven waves. Due to the jump in the wave veloc- 
ity at the interface, the wavenumber increases, i.e. the resolution 
requirements in the core are more restrictive. This holds, in partic- 
ular, for weak magnetic fields, where the jump in the wave velocity 
is larger, and for large / modes with higher frequencies. 
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The damping time at a given magnetic field strength B > 
10 13 G varies only by at most a factor of two between different 
I, but may vary significantly with the magnetic field strength. This 
indicates that the damping is dominated by the magnetic field and 
does not depend sensitively on the mode structure itself, i.e. there 
will always be a part of the continuum that is able to drain the en- 
ergy from the crustal oscillations. 

A more detailed discussion of the damping of the different 
crustal modes is beyond the scope of this work. Our results do not 
favor a crustal-mode interpretation of the observed QPOs in SGRs, 
because any crustal shear mode is damped sufficiently fast for mag- 
netic field strengths well below the typical magnetar field strengths 
~ 5 x 10 13 G. Although there might exist SGRs with weaker mag- 
netic fields (Rea et al. 2010 ), QPOs have only been observed so far 
in magnetar s with the strongest fields. 



5.2.2 n > shear modes 

The higher radial overtones (n > 0) of the shear modes have fre- 
quencies above 500 Hz (see Table[2| and are usually used to explain 
the QPOs of SGR 1806-20 with the highest frequencies of 625 and 
1840 Hz. As these modes have at least one node inside the crust 
computing their evolution demands much higher spatial resolution 
than necessary for the n — modes. Thus, it was practically im- 
possible for us to follow their evolution over several Alfven cross- 
ing times to the same accuracy as for the n — modes. However 
this would have been necessary to draw more reliable conclusions 
about damping times or interaction of these modes with the Alfven 
continuum in the core. Nevertheless, we can make some qualitative 
statements. 

The time needed for a shear wave (vs ~ 1000 km/s) to travel 
through the crust ( Ar cr ust ~ 1 km) corresponds to the inverse of 
the frequency of the first overtone (n — 1), which is of the or- 
der of / rsj 1 kHz. Assuming that the wave travels inside the crust 
along the ^-direction (travel path Ar ~ 107rkm), a similar esti- 
mate results in frequencies vs / Ar ~ 30 Hz which is of the order 
of the frequency of the fundamental n = oscillation. Therefore, 
we may conclude that the n — modes represent waves that travel 
predominantly parallel to the crust-core interface, while the n — 1 
modes correspond to waves that travel radially. This may explain 
the strong dependence of the n — modes on the angular number 
/ and the weak dependence of the n — 1 modes (see Table [2] with 
the corresponding discussion, and |Sotani et al.||2007| >. Hence, we 
expect the n — modes, derived with isotropic shear modulus, to 
be much more affected by the presence of an anisotropic magnetic 
field than the n = 1 modes. In particular, near the equator where 
the coupling between the crust and the core is weaker than close 
to the pole (see Eq.[42]> the shear waves may travel back and forth 
in the crust without interacting strongly. Additionally, at the equa- 
tor the magnetic field is almost parallel to the ^-direction, i.e. the 
direction of the Alfven waves is perpendicular to the direction of 
the n — 1 shear waves. This suggests that the n = 1 shear waves 
are not influenced strongly by the presence of moderate magnetic 
fields (< 10 15 G). According to these theoretical considerations we 
expect the overtones of the shear modes to survive longer than the 
n = modes. 

To investigate this issue, we have performed two simulations 
with our reference model at magnetic field strengths of 2 x 10 14 G 
and 5 x 10 14 G. The resolution for antisymmetric simulations (/ = 
{2, 4, 6, ..}) was 150 x 50 zones. The initial perturbation consisted 
of the n = 1, I = 2 mode of the crust only. In Figure [5] we show 
the main contributions to the oscillations as obtained by Fourier an- 
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Figure 5. The spatial structure of the first n = 1, antisymmetric crustal 
shear modes in the presence of a moderate magnetic field of 5 x 10 14 G. 
The oscillation patterns are strongly distorted from the typical £ -dependence 
of the spherical harmonics and only exist in regions close to the equator, 
where the coupling between the core and the crust is weakest. We therefore 
label their angular dependence with The dashed lines indicate the region 
of the crust and the color scale ranges from white (minimum) to red-black 
(maximum). 



alyzing the time evolution of the simulation with 5 x 10 14 G. The 
QPO patterns are compressed towards the equator and strongly dis- 
torted from the typical /-dependence of spherical harmonics of the 
purely shear eigenmodes of the crust (see Appendix [5]). To account 
for this difference we label these QPOs with /'. Naturally, the initial 
data of the undistorted /-mode excites many distorted /'-QPOs. As 
expected from the theoretical considerations above, the strongest 
amplitudes of the oscillations appear near the equator, where the 
coupling to the core is weakest (see Fig.[5]). The magnetic field also 
increases the spacing between the frequencies of successive n = 1, 
I' modes from Af « 1 Hz without field to Af « 10 Hz (see |Sotani 
|et al. 12007] for a discussion of purely shear eigenmodes) 



To study the behavior of the different QPOs, we calculate the 
corresponding overlap integrals (Eq. ( |B11| )) but taking the spatial 
structure obtained from the Fourier analysis (Fig. [5} as basis func- 
tions. The time evolution of the maxima of these overlap integrals 
for /' = {2,4,6,8} can be seen in the left panel of Figure[6] Indeed 
all of /'-QPOs are excited by the purely shear n = 1, / = 2 eigen- 
mode perturbation. Since the frequencies of the different n = 1 
modes are very similar, i.e. hard to disentangle in any analysis of 
the simulations, we average all n = 1 modes to estimate the to- 
tal damping time of the n = 1 QPOs. To this end we calculate 
the overlap integral with the radial function Hi (r, 6) — R\ r cor- 
responding to the pure shear eigenmodes obtained with Eq. ( |B8| >. 
This effectively averages over the angular dependence, and pro- 
vides a measure of how strong the ensemble of all /', n — 1 
QPOs is excited. The corresponding plots for B = 2 x 10 14 G 
and B — 5 x 10 14 Gare shown in the right panel of Fig. [6] As 
indicated by the fitting functions with damping times of 150 and 
230 ms, the damping timescale of the n — 1 QPOs is much longer 
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Figure 6. Left panel: The normalized overlap integrals of the evolution for the first few n = 1, antisymmetric magneto-elastic modes at B = 5 x 10 14 G. The 
approximate mode structure is taken from the Fourier analysis (see Fig. [5}. Strictly speaking there are no modes in terms of the result of the linear analysis, but 
the influence of the magnetic field is insufficient to destroy the coherent oscillations in all parts of the crust. The initial perturbation of the n = 1, 1 = 2 crustal 
mode excites a large number of the magneto-elastic modes which are all damped on a timescale of ~ 150 ms. Right panel: The overlap integral performed for 
a basis with purely radial dependence S^(r, 6) = R\ r (see Appendix|Bj of the n = 1 modes. These integrals are a measure of how strong the ensemble of 
all/ 7 , n = 1 modes is excited. The damping time, indicated by the fits to an exponential (dashed lines), for 2 x 10 14 G (5 x 10 14 G) of r = 230 ms (150 ms) 
is much longer than that of the n = modes, which is of the order of several ms only. 



than for the n — modes at the given magnetic field strength of a 
few 10 14 G. 

However, one has to be very cautious at this point. With the 
resolution used here, we are not able to resolve the Alfven oscil- 
lations inside the core of the neutron star which could damp the 
crustal shear oscillations resonantly. At 5 x 10 14 G the fundamen- 
tal Alfven oscillation is about 2 Hz. To resolve the resonant cou- 
pling to the crustal n — 1 mode of roughly 500 Hz one would need 
the 250th overtone. Simulations with appropriate grid resolution 
would take of the order of years. For the damping process itself this 
lack of resolution should not be a problem, because the numerical 
method employed should take all necessary information into ac- 
count, i.e. the Riemann solver at the crust core interface considers 
all the local information of possible waves traveling into the core. 
The problem arises inside the core, where a low resolution leads 
to an averaging out of all fine- scale structure. Hence, the energy of 
the Alfven overtones of the continuum is transformed to resolved 
low order oscillations, i.e. we cannot trust the oscillations inside the 
core. However, as long as the Alfven oscillations do not reach the 
crust at the opposite side of the star at £ « £a = 0.5/ f a = 1-2 s 
(for B = 5 x 10 14 G),a simulation of the crust region should give 
correct results. Therefore, the present estimate of the damping time 
of the n = 1 overtones should be considered as a lower limit, be- 
cause the main effect we are missing is the excitation of crustal 
magneto-elastic QPOs by incoming Alfven oscillations of the core. 

5.2. 3 Results for different equilibrium models 

All previous results were obtained for a model based on the 
APR+DH EoS and a mass of 1.4 M©. Simulations using other EoS 
or different masses for the equilibrium model yield qualitatively 
similar results. 

In table 2 of [Gabler et al."|( |2011| ) we showed that the damp- 
ing of purely crustal shear modes occurs on the order of 100 ms or 
less for a variety of EoS and masses at a magnetic field strength of 
B — 10 14 G. The initial models that we used in those simulations 
employed a different prescription for the internal energy (ideal gas) 
than the values provided in the EoS tables. As these models were 



EoS 


n = 0, 1 = 


r [ms] at B = 10 14 G 
2 n = 0, 1 = 3 n 


= 0, 1 = 9 


APR+DH 1.6 


45 


47 


21 


APR+DH 2.0 


38 


41 


17 


L+DH 1.6 


57 


61 


27 


L+DH 2.0 


52 


56 


25 



Table 4. Damping timescales r due to resonant absorption of crustal shear 
modes by the Alfven continuum for initial perturbation modes I = 2, 1 = 3, 
and I = 9 for different combinations of equations of state at B = 10 14 G. 
The number in the labeling of the EoS represents the mass of the neutron 
star model in Mq . 



thus thermodynamically inconsistent, we have recomputed them 
with the appropriate internal energy values. The resulting damp- 
ing times are at most 20% longer than the values reported in Gabler 
|et al.|(|2011| ). Some of the corrected damping times are shown in 
Table|4] The conclusions drawn in Gabler et al. (2011) are not af- 
fected at all by this change, because the maximal damping time is 
still less then 100 ms in all cases. 

The variation of the damping times with the EoS at a given 
magnetic field is not surprising. The relative size of the crust of 
these models varies roughly by a factor of 3 (see Sotani et al. 2007), 
and the shear modulus of both crustal EoS is of comparable size. 
This explains the smallness of the observed variations in Gabler 
|etal. 1(2011) which do not exceed a factor of 5. A significantly lower 
shear modulus, as proposed in Steiner & Watts ( 2009), would lead 
to even shorter damping times of the crustal shear modes. 

The influence of the details of the EoS are not substantial be- 
cause when trying to explain the frequencies of the QPOs observed 
in SGRs as shear oscillations, the shear modulus should lie in the 
range we use in this work. Otherwise it is already impossible to 
reproduce the correct range of frequencies within the crustal oscil- 
lation model. Changing the shear modulus somewhat would have 
only a modest effect on the damping times. However, even for a 
hypothetically exotic shear modulus, which could be one order of 
magnitude larger than the actual values we use, the crustal shear 
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Figure 7. The Fourier amplitude inside the neutron star for the model 
APR+DH 1.4 at B = 4 x 10 14 G. Shown are the first two lower QPOs 
l[, ±} and L^, the first four upper QPOs [/* (+) , , , and u[~ ] 



and some selected edge QPOs e[ + \ E^~\ and E^~ J . (The figure for 
E^^ was not very clear due to contamination with other QPOs.) The plus 
and minus sign indicate symmetry (+) and antisymmetry (-) with respect to 
the equatorial plane. Magenta lines indicate magnetic field lines, and black, 
dashed lines the location of the crust. The color scale ranges from white- 
blue (minimum) to red-black (maximum). 



*(+) 



oscillations would be damped much too fast to explain long-lived 
QPOs. We therefore argue that we can safely exclude shear oscil- 
lations as a viable explanation of observed magnetar QPOs, for the 
magnetic field configurations studied here. 



5.3 QPO structure for intermediate field strength 

In |Gabler et aT] ( [20T0l|20TTl > we have shown that the crust signifi- 
cantly changes the structure of the QPOs for intermediate magnetic 
field strengths between 5xl0 13 — 10 15 G compared to models with- 
out crust. In this section we investigate the behavior of the Alfven 
QPOs of the core in this regime further. 



5.3.1 QPO structure at 4 x 10 14 G 

In this subsection we will analyze the results of two simulations 
with a resolution of 100 x 40 zones and computational domain 
[0, r s ] X [0, 7r/2]. We perform one simulation with / = 2 initial data 
and antisymmetry with respect to the equatorial plane, and another 
symmetric one with / = 3 initial data. Both runs were evolved up 
to t ~ 5 s. 

For the model APR+DH 1 .4 with a magnetic field strength of 
B = 4 x 1Q 14 G we fin d three different families of QPOs. Follow- 
ing 



Cerda-Duran et al. 



(2009) we label the lower QPOs as lL ±} and 



the upper QPOs as Un • A new family of QPOs appears, which 
we call edge QPOs and label them as To avoid confusion 

with the previous work of Cerda-D uran et al.| ( |2009| ), and because 
the fundamental upper symmetric QPO has special properties as we 



will show below, it will be labeled as at low magnetic field 
strength. The plus and minus sign in the description of the QPOs 
indicate symmetry (+) or antisymmetry (-) of the QPO with respect 
to the equatorial plane. 

In Fig. [7] we plot the distribution of the Fourier amplitude in- 
side the star. Lower QPOs, as shown in panels a and b, are attached 
to field lines which close inside the core. In contrast, upper QPOs 
are located closer towards the magnetic poles at open field lines 
(panels c, d-f). A third family, the edge QPOs, are connected to the 
open field line inside the core of the neutron star which just fails to 
close inside (last open field line). These QPOs can be seen in pan- 
els g-i of Fig.[?] The two QPOs and E^ have very similar 
frequencies. Due to limited evolution time and hence limited reso- 
lution for the Fourier transform both QPOs contribute significantly 
to the Fourier signal at the corresponding frequency as can be seen 
in panel i. Similarly, the figure for E^ contains some contribution 
of Uq^ along the field lines crossing the equator between a radius 
of 1 and 4 km. In both panels the edge QPOs are concentrated on 
the field lines which cross the equatorial plane at around 5 km. The 
naming of the edge QPOs becomes clearer in Fig.|8] where we plot 
the Fourier amplitude of the velocity, averaged per field line in the 
frequency-radius plane. The maxima indicate the position of the 
QPOs. The red and green lines are the continuum of frequencies 
obtained with the semi- analytic model introduced in Cerda-Duran 
|et al.| ( [2009] > and adopted here to the problem in the presence of the 
crust. 

The different families of QPOs mentioned above and shown 
in Fig.[7]can also be identified in Fig.[8] The lower QPOs are at- 
tached to the closed field lines which cross the equatorial plane 
near 6.5 km. Since they are connected to the closed field lines these 
QPOs do not have a preferred symmetry, and hence are present 
in symmetric and antisymmetric simulations. The continuum of 
frequencies derived with the semi- analytic model (green lines in 
Fig.[8j has a minimum at the point where we find the lower QPOs . 



We thus interpret the as turning point QPOs (see 
With the exception of the fundamental ui + \ which is located al- 
most at the turning point of the semi-analytic model, the upper 
QPOs are localized in the continuum of the open field lines (red 
lines in Fig.[8j that cross the equatorial plane at 2 — 4 km. 

The members of the new family of QPOs, also obtained in 
|Colaiuda et al.| ( [2009] >, are called edge QPOs because of their posi- 
tion in Fig. [8] These QPOs are related to those parts of the contin- 
uum, obtained with the semi-analytic model (red lines), which do 
not connect to the continuum of the closed field lines (green lines) 
(see also the sketch in Fig.[TJ. For more details on the interpretation 
of the QPO structure without crust we refer to Cerda-Duran et al. 
( [2009} and with crust to |Gabler et ak] ( [20T0l ). 



Levin 



2007). 



5.3.2 Differences caused by the presence of the crust 

The lower QPOs attached to the closed field lines are reproduced 
qualitatively similar as in the case without crust. The only differ- 
ence is that they are limited to the field lines which close inside the 
core and do not extend into the crust. 

However, the upper QPOs which are located near the pole 
in models without crust ( Cerda-Duran et al. 2009, Colaiuda et al. 
2009 ; Sotani et al. 2007]), can now be found at substantial distance 
from the poles, i.e. at lower latitudes (see also figure 4 in Gabler 
|et al.|20lT) . In simulations without crust the oscillations were as- 
sociated with the maximum at the turning point of the continuum 
at the pole, while if the crust is included, we obtain the maximal 
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Figure 8. The averaged Fourier amplitude along different magnetic field 
lines labeled by the radius where they cross the equator for model APR+DH 
with M = 1.4M© and B = 4 x 10 14 G. Red and green lines give the 
Alfven continuum obtained with the semi-analytic model. The locations of 
the QPOs are indicated by black dots. Left panel: antisymmetric simulation; 
right panel: symmetric simulation. The color scale ranges from white-blue 
(minimum) to orange-red (maximum). 



amplitudes away from the pole and inside the continuum predicted 
by the semi-analytic model in the absence of a crust. One possible 
interpretation of this new feature is that the shear modulus in the 
crust alters the propagation of magneto-elastic oscillations in the 
region near the pole in such a way that standing waves cannot form 
at all along individual field lines or, if they form, they go quickly 
out of phase with nearby field lines in this region. 

At this point it is helpful to recall the problem of the reflection 
of plane -parallel waves, where the reflection coefficient depends on 
the jump in the propagation velocity or equivalently on the index of 
refraction. The stronger the jump in the index is, the larger is the 
fraction of the incident wave which becomes reflected. Using this 
analogy we would expect that the smaller the difference in prop- 
agation velocity at the crust-core interface is, the more refraction 
into the crust should occur and less reflection back into the core 
should be produced. If significant refraction into the crust occurs, 
no stable standing waves can be maintained during the evolution. 
When following the crust-core interface from the pole towards the 
equator the coupling between crust and core becomes weaker, see 
Eq. ( |42| where the dependence of the coupling factor is realized in 
b r . The magnetic field in the radial direction, and hence the Alfven 
velocity, decreases with increasing 6. Therefore, the jump in the 
propagation velocity increases, and thus the fraction of the wave 
which is reflected will also increase along this trajectory. For suffi- 
ciently low magnetic fields there should always be a region near the 
equator were almost perfect reflection occurs. However, when fol- 
lowing the crust-core interface from the equator towards the pole, 
one will reach a characteristic magnetic field strength where insuf- 
ficient reflection occurs to maintain stable standing waves along a 
certain field line. At this point we find the maximum amplitude of 
the QPO in the presence of a crust. We will investigate the reason 
for this behavior further in Sec. 15.3.41 

The new position of the maximum amplitude of the QPOs is 
thus determined by two effects. Near the pole the magnetic field 
lines get out of phase due to the interaction through the extended 
crust, because a significant fraction of the oscillation is refracted. 
The magnetic field lines can be seen like strings which are not at- 
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Figure 9. Evolution of the magnetic plus kinetic energy per field line di- 
vided by the sum of the energy over all field lines. The field lines are la- 
belled by their crossing point with the equator. Note that as the total energy 
decreases with time the apparent increase of the energy of the lower QPOs 
is due to color rescaling at every time step. Their amplitude actually de- 
creases, but more slowly than that of the upper QPOs. The left panel shows 
antisymmetric and the right panel symmetric simulations. The color scale 
ranges from white-blue (minimum) to orange-red (maximum). 



tached to a rigid but rather "moving" boundary, i.e. the crust which 
responds to the oscillations. The resulting effect is a strong cou- 
pling of different magnetic field lines and, consequently, an energy 
transfer between the lines due to the scalar shear modulus. In this 
region and for magnetic field strengths studied here, each field line 
seems to act as a damped oscillator. Near the equator the magnetic 
field lines get out of phase due to phase mixing like in the case 
without extended crust. The additional damping close to the pole 
makes the corresponding upper QPOs to be shorter lived than in 
the case without crust. 

This effect can be observed in the right panel of Fig.[9] where 
we show the magnetic plus the kinetic energy per field line divided 
by the total magnetic plus kinetic energy at the given time as a 
function of time for different field lines. The initially excited QPOs 
attached to the field lines between 4.5 and 5 km disappear rapidly 
after about 0.5 s. In contrast, the lower QPOs and the fundamental 
symmetric QPO near the pole persist during the whole evolution. 
Fig.[9]may suggest that the lower QPO gain energy with time. How- 
ever, this apparent energy increase is not a physical effect, because 
the total energy decreases with time due to numerical dissipation. 
Hence, the relative amplitudes of the energy of the lower QPOs in- 
crease, while their absolute amplitude decrease slightly because of 
numerical dissipation. 

Compared to simulations without crust ([Cerda-Duran et aL] 



2009 



we find a new fundamental upper QPO U* +) . This QPO ap- 
pears because the boundary condition at the crust-core interface 
causes a reflection which results in a node at this surface. With- 
out crust the boundary condition at the surface of the star im- 
plies a maximum there and the fundamental oscillation has the 
node at the equator in this case. Therefore, the symmetric QPO 
^o + — must nave an additional node inside the core (see panel e in 

r(+) 



10 1. £7* situated between 1 and 



Fig.^or the second row of Fig 

2 km (Fig.|9] the right panel) decays less rapid than the other up 
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per QPOs, \ U^ , and £/£ >0 . This may be related to the fact, 

that at B = 4 x 10 15 G, is located close to the maximum of 
the continuum (see Fig.[8j. There, the gradient of the continuum is 
less steep, and neighboring field lines get out of phase less rapidly. 
This behavior is similar to a turning point QPO, which persists for 
longer time than edge QPOs ( Levin 2007| ). 



r(-) 



5.3.3 Changing QPO position with increasing magnetic field 

For magnetic field strengths between 10 14 and 10 15 G the simu- 
lations reveal that the location of the upper QPO, changes 
within the neutron star. This was not observed in the case of pure 
Alfven oscillations in [Cerda-Duran et al.] ( |2009| > or |Sotani et al.| 
(008b ), where the upper QPOs were always observed close to the 
pole. Fig.[l0] (first three columns) shows the new effect, where we 
plot the spatial structure of the Fourier amplitude of the QPOs. 
When increasing the magnetic field from 2 x 10 14 to 10 15 G the 
upper QPOs move from a location near the pole towards the 



equator. The change in position of the is shown in Fig 

only for and u[~\ but holds for all higher overtones as well 
and does not depend on the symmetry. For the fundamental U*^ 
the dislocation is less (upper row in Fig.[l0j. One can understand 
this behavior at least partially with the help of the semi-analytic 
model. Fig.|8] shows that the frequencies and the symmetry of the 
QPOs are correctly predicted, if we assume that the Alfven wave 
is reflected at the crust-core boundary. However, the semi-analytic 
model cannot explain where within the continuum the QPOs are 
situated. 

Remembering the analogy with the reflection of plane-parallel 
waves in the preceding subsection and bearing in mind that with 
increasing magnetic field b r the relative jump in the propagation 
velocity on both sides of the crust-core interface at a given position 
decreases, more parts of an incident wave get refracted into the 
crust. This means that the point where stable standing waves can 
be maintained should move towards the equator. This is exactly 
where for magnetic fields < 10 15 G the QPOs 



10 



10 



confirmed in Fig. 
move from close to the pole towards the equator, as the magnetic 
field strength increases. 



5.3.4 Reflection of pulses and spread in crust 

In the absence of an elastic crust, Alfven wave packets are supposed 
to travel approximately along magnetic field lines, as the character- 
istic direction of propagation of any magnetic perturbation coin- 
cides with the direction of the magnetic field. However, when a 
crust is added this picture changes. The direction of propagation 
of magneto-elastic waves no longer coincides with the magnetic 
field direction (compare the different eigenvalues in this case given 
in Eq. p7]>). One would therefore expect a perturbation, traveling 
along magnetic field lines from the center of the star towards the 
surface to spread out past the crust-core interface. Such a spread 
is strong for low magnetic field strengths, when the isotropic shear 
modulus dominates in the crust region and weak for high magnetic 
field strengths, when the opposite is true. This behavior is shown in 
Fig.JTT] where we display the renormalized sum of the kinetic and 
magnetic energy per field line for simulations at 5 x 10 14 G. The 
initial perturbation is restricted to a limited region of the star about 
4 km above the equator. The left panels show the expected behav- 
ior for a pulse which travels along a field line with no crust present. 
Two initial perturbations differing only by the location of the star 
get reflected at the surface and travel back towards the center. Any 




equatorial radius [km] 



equatorial radius [km] 



Figure 11. Evolution of the magnetic plus kinetic energy per field line di- 
vided by the sum of the energy of all field lines for the model APR+DH 
1.4. The field lines are labeled by their crossing point with the equator. The 
pulse reaches the surface at around 70 ms for the first time and the crossing 
time is about 130 ms. Left panels: simulation without crust, right panels: 
simulation including crust. The upper and low panels differ by the loca- 
tion of the initial perturbation. Color scale ranges from white (minimum) to 
red-black (maximum). 



deviations from traveling perfectly along the initially excited field 
lines is caused by a numerical coupling of different field lines and 
the very weak coupling through the boundary condition at the sur- 
face of the star. Taking a crust into account but imposing the same 
initial perturbations the wave packets are spread whenever entering 
the crust, which happens at around 70 ms for the first time, and sub- 
sequently after about every 130 ms (right panels). For initial data 
located at field lines crossing the equator around 2 to 3 km (lower 
right panel), the spread is more drastic. After some reflections there 
are phases (around 700 and 850 ms), when no significant perturba- 
tion amplitudes can be found around the field lines which initially 
carried the perturbation. 

However, the scenario of an initially localized wave packet 
considered here cannot rule out the existence of standing waves 
along individual field lines. Nevertheless, it suggests that additional 
effects may be introduced by the spreading of wave packets in the 
crust which probably change the Alfven continuum of our semi- 
analytic model, where one assumes standing waves which get re- 
flected at the crust-core interface (or at the surface for stronger 
magnetic fields). 

5. 3. 5 Conservation of angular momentum 

Analyzing the convergence properties of our numerical simulations 
in the intermediate magnetic field case is more complicated because 
the contributions to the stress-energy tensor from the magnetic field 
and the shear are of the same order of magnitude. The two ex- 
treme regimes have been tested above (purely shear oscillations) 
or in Cerda-Duran et al. (2009) (purely Alfven oscillations). In the 
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2xl0 14 G 6xl0 14 G 10 15 G 1.5xl0 15 G 2 x 10 15 G 3 x 10 15 G 5 x 10 15 G 8 x 10 15 G without crust 




Figure 10. Structure of upper QPOs at different magnetic field strengths for the model APR+DH 1.4. The panels show the fundamental symmetric U* 
(upper panels), the symmetric (middle panels) and antisymmetric u[~^ QPOs (lower panels). The frequencies of the QPOs shown here are given in 

Table[5] The color scale ranges from white-blue (minimum) to red-black (maximum). 

present approach, the total angular momentum J to t = f S^dV 
is the only globally conserved quantity, with dV = ^drdOdcj). 
We do not expect conservation of the energy of the perturbation in 
our simulations, because neglecting the coupling to poloidal oscil- 
lations and assuming purely poloidal magnetic fields, renders the 
deviations of the total energy from the energy of the unperturbed 
background configuration to be of second order in the perturba- 
tions, while our approach is accurate to first order. 

The total angular momentum J to t should be conserved in- 
side the computational volume, but the boundary condition we 
have chosen (see Sec. |2.3| > allows for non-vanishing flux through 
the surface. As we chose initial perturbations with the angular de- 
pendence of the vector spherical harmonics the angular momenta 
in both hemispheres cancel by construction and the total angular 
momentum of the star is zero. For antisymmetric perturbations the 
losses/gains through the surface cancel respectively, while for sym- 
metric perturbations there remains a non-zero contribution. 

Fig. [12] shows the variation of the total angular momentum 
during the evolution for an symmetric simulation with / = 3 ini- 
tial data at three different grid resolutions. When analyzing the dif- 
ferences between the curves (see Section [5T| we obtain the order 
of convergence of 1.95, which is near the expected second-order 
convergence. To estimate the absolute magnitude of the resulting 
angular momentum error we compare our perturbation to a rigidly 
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rotating sphere with the same total angular momentum. Taking the 
typical total angular momentum during the simulation, and com- 
paring it to that of a rigidly rotating sphere J — 2/5MR 2 Q — 
2/5MRv TOt , where Q = v TOt /R is the rotation frequency, we ob- 
tain a maximal velocity which is only a fraction of the perturbation 
used in the simulations v ro t/v pe rt ~ 10 -10 . Thus, the total angu- 
lar momentum introduced by our perturbation is very small. More- 
over, the total angular momentum variations converge to zero, i.e. 
compared to the numerical errors the losses of angular momentum 
through the surface is a small effect not affecting our simulations. 



5.4 QPO structure for high magnetic field strength 

For strong magnetic fields we are interested in the structure of 
magneto-elastic QPOs, but not in the damping of crustal modes, 
as the latter are damped already at much lower field strengths (see 
Sec. |5.2| ). Therefore, we can reduce the grid resolution, i.e. the com- 
putational costs. In the strong field case we thus use a uniform ra- 
dial grid with 100 zones and an angular grid with 80 zones in the 
interval [0, it]. 

For very strong magnetic fields, B > 5 x 10 15 G, the maxima 
of the Fourier transform align towards the polar axis. With increas- 
ing magnetic field the influence of the shear inside the crust be- 
comes negligible, and the QPO pattern approaches that expected for 
the purely magnetic limit (see figure 3 in Cerda-Duran et al. 2009), 
in agreement with the semi-analytic model. Another effect caused 
by the anisotropy of the shear modulus is a more wide- spread spa- 
tial structure of the QPOs compared to the case without crust. This 
can be inferred from the last two columns of Fig.[l0] where QPOs 
are still quite extended inside the crust (B = 8 x 10 15 G), and in 
the model without crust. 

Between the two extremes, the QPOs are confined in the core 
(B < 10 15 G). For strong magnetic fields (B > 10 15 G) there is 
a transition from the QPO structure observed in Section [53] to the 
purely magnetic case (see Fig. [To| from the 3 - 6 column). Between 
10 15 and 2 x 10 15 G the QPOs begin to have significant amplitudes 
in large parts of the crust. 

This transition becomes clearer in Fig. [13] where we plot the 
Fourier amplitude for individual field lines averaged over the length 
of the line and labelled by their crossing point with the equatorial 
plane. The solid lines represent the continuum as obtained by the 
semi-analytic model, where we assume reflection at the crust-core 
interface in the upper row and reflection at the surface of the star 
in the lower row. The obtained frequencies are very similar in both 
cases, because the travel time of the waves is dominated by the time 
spent in the core. We note a change of the structure of the QPOs 
with the different boundary conditions. We have already noted in 
Fig.[l0|that in the case of weak magnetic fields the QPOs move 
from being near the pole towards the equator, for increasing mag- 
netic field strength below 10 15 G. The same can be observed in the 

. Between 10 15 
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first two panels of Fig. 

and 2 x 10 15 Ga transition occurs from reflection at the crust-core 
interface (for weaker magnetic fields) and reflection at the surface 
of the star (for stronger fields). Therefore, we do not expect any 
of the two approximations to agree with the semi-analytic model. 
Neither can we assume that the oscillations get reflected only at the 
crust-core boundary nor only at the stellar surface. When increas- 
ing the magnetic field beyond 2 x 10 15 G the maximum amplitude 
of the QPO aligns again towards the polar axis, and the numeri- 
cally obtained pattern approaches the expected structure, which is 
similar to the case without the crust. The main differences are the 



presence of the fundamental U*^ at finite frequency and a broader 
maximum because of the anisotropic shear contribution. The QPOs 
are moving from the pole towards the equator according to the fre- 
quency predicted by the semi- analytic model for reflection at the 
crust-core interface, but when the effects of the magnetic field be- 
come comparable to those of the shear modulus in the crust, the 
QPOs reach the end of the continuum and jump from the symmetric 
(antisymmetric) branch to the antisymmetric (symmetric) branch of 
the next part of the continuum. 

What happens to U* + \ which has no possibility to jump to? 

For strong magnetic fields, B > 2 x 10 15 G has different 

features. First, there is no node along the field lines, due to the 

change of the boundary conditions, which require a maximum at 

the surface. Second, as can be seen in the last few columns of the 

upper row of Fig.[l0] there are nodes perpendicular to the field lines 

at a given frequency, and we even find two different contributions 

to the at slightly different frequencies. Between 2 x 10 15 G 

and 5 x 10 15 G there is one node perpendicular to the field lines, 

while for B ^ 6 x 10 15 G the QPO has predominantly two nodes 

in that direction. Note that both features are always present for 

B > 2 x 10 15 G. This splitting can also be seen in the leftmost 

three panels of the lower row in Fig.[l3] The corresponding fea- 

— r\ 
tures are located at frequencies below the fundamental, Uq . The 

same panels also show, that with increasing magnetic field strength, 

the relative Fourier amplitudes and the frequencies of those features 

decrease with respect to the fundamental frequency. 

The picture of the transition from one asymptotic behavior (re- 
flection at the crust-core interface) to the other (reflection at the 
surface) gets supported by the frequencies obtained for the differ- 
ent QPOs shown in Table[5] There we compare the frequencies cor- 
responding to the maxima of the Fourier amplitude with those of 
the fundamental oscillation obtained with the semi- analytic model. 
Thus, we use the version with reflection at the crust-core interface 
up to 10 15 G, while we set the boundary at the surface of the star 
for stronger fields. For magnetic fields up to B < 8 x 10 14 G, 
U*^ has a similar frequency as the fundamental obtained with the 
semi-analytic model. The frequencies of the other QPOs approxi- 
mately behave as 1 (£/* (+) ) : 3 (£/ (+) ) : 5 (U[ +) ) and 2 ([/ (_) ) : 
4 (U[~ } ) : 6 (t/ 2 (_) ). For stronger magnetic fields B > 3 x 10 15 G 
the frequencies ratios approach 1 (Uq~^) : 3(11^) : 5 (U^) 
and 2 (U^ 



4([/-[ + ' ) ). The two asymptotic integer relations 
between successive overtones and their order is what is expected 
from the semi- analytic model in the two regimes. In the intermedi- 
ate regime 8 x 10 14 < B < 3 x 10 15 Gthe frequencies change 
smoothly from one relation to the other. 

Deviations from exact integer ratios may have different rea- 
sons. First, the time of numerical integration is limited. Therefore, 
the spectral resolution of the Fourier analysis is limited, too. Sec- 
ond, for the lowest magnetic field shown here, B ~ 2 x 10 14 G, 
not all upper QPOs have reached their position near the polar axis, 
i.e. their frequencies still lie in the continuum, resulting in lower 
frequencies. Third, in particular in the transition regime it is some- 
times difficult to identify where the maximum of a QPO is located. 
The interesting QPO may be excited only very weakly by our ini- 
tial data, and/or some other QPO may be excited more strongly at 
a similar frequency. This occurs more frequently for higher over- 
tones, because there the different continua overlap (see Fig.[8]or 

ED- 
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Figure 13. The averaged Fourier amplitude along different field lines labelled by the radius where they cross the equatorial plane for different magnetic field 
strengths. In the last panel we show the expected behavior for a model without crust. The y-axis represents the frequency rescaled to the fundamental QPO 
frequency predicted by the semi-analytical model (see Table|5j. The color scale ranges from white-blue (minimum) to orange-red (maximum). Black and green 
lines represent the continuum when taking the crust into account. Dashed, red lines indicate regions where the semi-analytic model is not supposed to work. 
We do not know where the waves get reflected and thus cannot identify whether to consider the lines as open or closed. 
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2 x 10 14 


1.0 


1.0 


1.9 


2.9 


3.6 


5.0 


5.7 


±0.24 


4 x 10 14 


2.0 


1.0 


1.8 


2.8 


3.5 


4.9 


5.5 


±0.20 


6 x 10 14 


3.0 


1.0 


1.8 


2.8 


3.1 


4.7 


5.0 


±0.10 


8 x 10 14 


4.0 


1.0 


1.8 


2.5 


3.5 


4.3 


5.1 


±0.20 


10 15 


5.0 


0.9 


1.3 


2.4 


3.3 


4.2 


5.3 


±0.10 


1.5 x 10 15 


7.4 


0.8 


1.1 


2.2 


3.2 


4.1 


5.1 


±0.10 


2 x 10 15 


9.8 


0.7 


1.0 


2.0 


3.0 


4.0 


5.1 


±0.10 


3 x 10 15 


14.7 


0.5 


0.9 


2.0 


3.0 


4.0 


5.1 


±0.10 


5 x 10 15 


24.4 


0.5 


1.0 


2.0 


3.0 


4.0 


5.1 


±0.10 


8 x 10 15 


39.1 


0.3 


0.9 


2.0 


3.0 


4.1 


5.1 


±0.10 



Table 5. The relation of the frequencies of the lowest QPOs for the model APR+DH 1.4, obtained by analyzing the local maxima of the Fourier amplitudes, to 
the fundamental frequency Uo, obtained by the semi-analytic method. Note that the fundamental of the semi-analytic model changes symmetry at a magnetic 
field of about 10 15 G, such that Uq ~ for B ^ 10 15 G and Uq ~ for B > 10 15 G. The last column shows the uncertainty in the separation 

of two successive frequencies in the Fourier spectrum. It gives an estimate of the error caused by choosing the position of the local maxima, but it does not 
include other numerical errors. 



5.5 Crustal modes in the gaps of the Alfven continuum? 

|van Hoven & Levin| ffiTT) and |Colaiuda & Kokkotas| ( [20TT] ) 
have pointed out the possibility of crustal modes, which may have 
frequencies outside of the continuum of the core. These modes 
would be only very weakly coupled to core oscillations, because no 



Alfven wave of any field line could match the necessary frequency. 
In the models we have studied here, we find gaps in the contin- 
uum only between the lowest overtones of Alfven oscillations (see 
Fig.[i}, e.g. for model APR+DH 1.4 already the continua of u[~^ 
and overlap and there are only gaps between Jj[ ~^ and U$ + \ 
Uq^ and Uq ~^ and below Uq~\ For the other models shown in 
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Figure 14. Frequency continua of model L+DH 1 .4 at magnetic fields be- 
tween B = 10 15 G and B = 4 x 10 15 G. Black, solid lines indicate 
the edges of the continua of the open field lines as indicated by the names 
and The horizontal red line gives the frequency of the / = 2 

crustal shear mode and the dashed red lines gives the frequency obtained 
by simulations with a free slipping crust (see text). The green shaded areas 
show continua having the same symmetry with respect to the equator as 
the / = 2 shear mode. The grey area shows the continuum with opposite 
symmetry. Magenta (blue) dashed lines show the range of magnetic field 
strengths where the frequency of the n = 0, I = 2 crustal mode lies in 
the gap between allowed Alfven continua for purely shear modes (for free 
slipping crust with magnetic field). 



the Table [2] the number of gaps is limited to a maximum of about 
three for models with the crust EoS NV and to two for models with 
the DH EoS. Note that we only consider the continua of the open 
field lines indicated by the black lines in Fig.[8] because they are 
decoupled from the continua related to the closed field lines (green 
lines). To have the fundamental n = 0, I = 2 oscillation of the 
crust in one of the gaps, very strong magnetic fields B > 10 15 G 
are required. 

As an example we choose the EoS L+DH 1 .4, because the cor- 
responding model has broader gaps than the APR+DH EoS. Fur- 
thermore, models with lower mass have shorter Alfven crossing 
times, and therefore weaker magnetic fields are necessary to have 
the shear mode in the gap between the lowest overtones of the con- 
tinua (see Table[7]). The spectral structure of this model is displayed 
in Fig.[l4] where we show the edges of the continua of the open 
field lines as a function of the magnetic field strength. The shaded 
areas between two edges represent the corresponding continuum, 
where crustal modes can be absorbed resonantly. The red line indi- 
cates the frequency 21.6 Hz of the purely shear, n = 0, 1 = 2 mode 
of the crust. Using model L+DH 1.4, we performed simulations 
with initial data with the crustal n = 0, I = 2 mode, allowing only 
for antisymmetry with respect to the equatorial plane. The oscilla- 
tions of the continuum associated with symmetric QPOs, as for ex- 
ample [7 (+) , are thus not allowed and cannot be excited. These for- 
bidden oscillations are indicated by the grey shaded region in Fig. 
[T4] We also confirmed that without imposing equatorial symmetry 
only Alfven QPOs having the same symmetry as the corresponding 
crustal mode can be excited with significant amplitude during the 
evolution. Therefore, the antisymmetric n = 0, 1 = 2 shear mode 
of the crust lies in the gap between and e[~^ for magnetic 
fields between 1.74 and 3.7 x 10 15 G (Fig.Jul Following previ- 
ous works by I van Hoven & Levin (2011) or Colaiuda & Kokkotas 
P011| > the crustal mode in the gap should not be damped signif- 
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Figure 15. Evolution of the overlap integral of the n = 0, / = 2 mode of 
the crust at different magnetic field strengths: B ={1.0, 2.0, 2.5, 3.0, 3.5, 
4.0} xlO 15 G. The dashed magenta line shows the expected oscillation of 
the purely shear I = 2 mode of the crust. All other lines show strong damp- 
ing of the excited I = 2 mode at early times (t < 20 ms). Contributions to 
the overlap integral at later times t > 20 ms originate from magneto-elastic 
oscillations. 



icantly, because there is no oscillation at the resonant frequency 
available in the continuum. 

However, we still find very strong damping of this crustal 
mode in our simulations, as indicated by the overlap integral of the 
n = 0, I = 2 shear mode in Fig. [15] Clearly the mode is damped 
after a few msec for all considered magnetic field strengths. Later 
contributions originate from coupled magneto-elastic pulses trav- 
elling through the whole star. Therefore, the time when they con- 
tribute to the overlap integrals depends inversely on the magnetic 
field strength, a behavior which is completely different from that 
of discrete oscillations for purely crustal shear modes. These ob- 
servations indicate that we do not have a weakly coupled system 
of two sub-systems (the crust and the core) but we are dealing with 
coupled, global magneto-elastic oscillations. To compare with what 
one would expect for an undamped purely crustal mode, we plotted 
the magenta, dashed line in Fig.[T5| 

We checked the influence of the magnetic field on the fre- 
quency of the purely shear mode by performing a series of simula- 
tions, where we apply an artificial boundary condition at the crust- 
core interface, i.e. we use the same condition = that would 
apply in the absence of the magnetic field (see also |Sotani et al.| 
2007). This allows the crust to slip freely on top of the core, and the 
oscillations inside the crust cannot be damped into the core. The 
frequency of the corresponding / = 2 crustal mode is influenced 
as displayed in Fig.[l4] This is expected, because global, poloidal 
magnetic fields of the order of B ~ 10 15 G and stronger begin to 
have measurable effects (Duncan 1998, Messios et al. 2001 ; Sotani 
|eTaLl|2007l |008al |Shaisultanov~ Eichler 2009}. Intuitively, one 
would expect the shear mode frequencies to increase with the mag- 
netic field strength, because the magnetic tension could be inter- 
preted as effectively augmenting the shear modulus (Messios et al.| 
|200T| |Sotani et al7| [2Q07 ). In contrast to |Sotani et al.| ( |2007| ) we 
find a decrease of the I = 2 mode frequency with increasing dipo- 
lar magnetic field strength (see Fig.[l4]). In their work |Sotani et al.| 
( 2007) neglected couplings between /- and / ± 2-modes due to the 
magnetic field, which led to discrete modes. However, by analyz- 
ing the evolution of the corresponding overlap integrals we observe 
strong excitation of the / = 4 mode by the / = 2 mode in our sim- 
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ulations. This strong coupling may explain the opposit change of 
frequency than expected from the study by Sotani et al. ( 2007). For 
a different magnetic field configuration also Messios et al. (2001) 
and |Sotani et"aL] ( |2007| > find that the frequency of the fundamen- 
tal/ = 2 crustal mode decreases with increasing magnetic field 
strength. 

If the frequency of the / = 2 crustal mode changes according 
to our simulations with the free slipping crust, the magnetic field 
strengths for which the mode lies in the continuum gap is limited 
to 1.4 x 10 15 to 2.2 x 10 15 G. In this regime we have performed four 
simulations at 1.5 x 10 15 , 1.7 x 10 15 , 2.0 x 10 15 , and 2.2 x 10 15 G 
(triangles in Fig.[l4]). In all of them we find the strong damping of 
the crustal mode. 

A second problem of matching of crustal frequencies into 
the gaps of the Alfven continuum emerges at the field strengths 
(> 10 15 G) at which we find the crustal frequencies in the contin- 
uum gaps, as there is no clear way of how to compute the Alfven 
frequencies. In this transition regime, reflection neither occurs pre- 
dominantly at the crust-core interface nor at the surface, there- 
fore, we do not claim that the continuum shown in Fig. [14] is per- 
fectly valid at all magnetic field strengths (compare the panels for 
1.5 x 10 15 G to those of 3 x 10 15 G in Fig.[l3}. However, by per- 
forming simulations for 10 different magnetic field strengths be- 
tween 10 15 and 4 x 10 15 G for the current model L+DH 1.4, we 
can ensure that the / = 2 crustal frequency lies in the gap be- 
tween the Alfven continua at least for one of the models. In none 
of the above simulations we find a different behavior than the one 
reported, i.e. we do not observe any of the crustal shear modes. 
For crust models with lower shear modulus, we expect the transi- 
tion to occur at lower magnetic field strength. The magnetic field 
required to fit a crustal shear mode into the gap between succes- 
sive continua decreases with decreasing shear modulus. However, 
a smaller shear modulus also lowers the relative importance of the 
shear terms compared to the magnetic field. 

In the example shown above, we may have just missed to 
match the frequency of the crustal mode to a gap of the continuum, 
as the frequency may have been changed due to the presence of the 
strong magnetic field, and because the continuum is probably not 
reliably predicted by the semi-analytic model. To this end we per- 
formed a large number (> 50) of simulations at different magnetic 
field strengths > 10 15 G and different equilibrium models but we 
never found any crustal shear mode at such strong magnetic fields. 

Generalizing the dipolar magnetic field configuration, which 
is our main model simplification in the current context, would prob- 
ably increase the complexity of the continuum, making it even 
harder to find gaps. However, there might arise new effects due 
to an entanglement of the magnetic field (see van Hove n & Levin| 
|20TT) . 

Crustal modes in the gaps of the Alfven continua of the core 
have been reported by Colaiuda & Kokkotas ( 201 1), but our results 
suggest that these QPOs are rather oscillations of the continuum. It 
is possible that the mode recycling technique (using the amplitude 
of the FFT at the frequency of the crustal shear modes) employed 
by |Colaiuda & Kokkotas |(|2011|) may give an e xcitation of the con- 
tinuum at this frequency. |van Hoven & Levin] ( |2012| > also observe 
gap modes. In their approach several 1 -dimensional eigenmodes of 
the magnetic field in the core are coupled to 2-dimensional eigen- 
modes in the crust. However, they do not prove that their assump- 
tion that the crustal dynamics can be described with an eigenvalue 
problem is valid also in the presence of a magnetic field. Thus, we 
expect that the ansatz of | van Hoven & L evin ( 2012) is valid only 
approximately for weak magnetic fields (B < 10 15 G). Such an 
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Figure 16. Normalized amplitude of the Fourier transform at the surface of 
the star (between pole and equator) as a function of the frequency for the 
stellar model APR+DH with 1.4 solar masses. The horizontal lines indicate 
the frequency for the upper QPOs predicted by the semi-analytic model. The 
color scale ranges from white-blue (minimum) to orange-red (maximum). 



eigenmode approach cannot describe the complicated coupling be- 
havior, including e.g. travelling wave packets, which would require 
the inclusion of a much larger number of coupled oscillators than 
those considered by | van Hoven & Lev in ( 2012). 



5.6 Threshold for the outbreak of the QPOs through the 
crust 



In Section [54| we noticed that for weak magnetic fields, B < 
10 15 G, the QPOs are largely confined to the fluid core, and that 
there exists a threshold magnetic field strength beyond which QPOs 
can be observed with significant amplitudes at the surface of the 
star. To quantify when magneto-elastic QPOs have a significant am- 
plitude in the crust of the neutron star, we studied their maximum 
amplitude at the surface. 

In Fig.[l6]we plot the amplitude of the Fourier transform at 
the surface for different magnetic field strengths as a function of 
the polar angle and the frequency. The color scaling is rescaled in 
each panel to enhance the main contributions at each magnetic field 
strength. We consider APR+DH 1.4 as our reference model. First 
we note, that for very strong fields (> 5 x 10 15 G), the frequen- 
cies approach those predicted by the semi-analytic model, because 
the influence of the crust should decrease for increasing magnetic 
field strength. There is some low frequency oscillation which cor- 
responds to U*^ (compare with Fig. 13), and there are some ad- 
ditional, strong Fourier modes near to 7r/2 resulting from the edge 
modes, which are stronger than in the case without crust. 

However, when decreasing the strength of the magnetic field 
from 5 x 10 15 G to 2 x 10 15 G, the Fourier amplitude of the QPOs at 
the surface decreases, e.g. the maximum decreases from 1.0 to 0.37 
in units normalized to the maximum amplitude at.B = 5xl0 15 G. 
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Model 


FFT(10 15 G) 
FFT (no crust) 


Model 


FFT (10 15 G) 
FFT (no crust) 


APR+DH 1.4 


0.07 


APR+NV 1.4 


0.06 


APR+DH 1.8 


0.11 


APR+NV 1.8 


0.005 


APR+DH 2.2 


0.33 


APR+NV 2.2 


0.04 


L+DH 1.4 


0.16 


L+NV 1.4 


0.016 


L+DH 1.8 


0.36 


L+NV 1.8 


0.03 


L+DH 2.2 


0.15 


L+NV 2.2 


0.009 



Table 6. The maximal magnitude of the Fourier transform at the surface 
of the star for a dipole magnetic field strength of 10 15 G. Simulations for 
different EoS, and with and without crust are considered and compared. 

As we have already seen in Fig.[l0| this is expected, because for 
decreasing field strength the crust will shield the QPOs. We fur- 
ther see that the additional structure in the Fourier amplitude at 
about the fundamental frequency of Uq~^ and the strong feature 
just below the frequency of close to the equator are increas- 
ing in amplitude relative to the upper QPOs. We note here that in 
the transition region of magnetic field strengths the correspondence 
between the frequencies of the semi-analytic model and the simu- 
lated QPOs should be taken with caution and we do not expect a 
simple QPO structure as in the two limiting cases. 

When decreasing the magnetic field strength towards 8 x 
10 14 G a dominant feature originates from Uq + \ In a sequence 
of similar plots for different magnetic fields not shown here, one 
can follow the slight change in frequency until reaching approxi- 
mately the value predicted by the semi- analytic model for reflection 
at the crust-core interface. Comparing with Fig.[l0|the correspond- 
ing spatial structure of the mode has still some amplitude inside the 
crust, i.e. our interpretation in terms of the QPO makes sense. 
The general trend is that the magneto-elastic upper QPOs exhibit a 
decreasing amplitude near the surface for decreasing magnetic field 
strength. 

For the weakest magnetic field (4 x 10 14 G), the upper QPOs 
have no strong amplitudes, because they are confined to the fluid 
core. The dominant contribution at the surface results from the edge 
modes in this case (see Fig.[7]). The QPO with the largest amplitude 
inside the crust is the edge QPO E^ (see Fig.[l6]) with a frequency 
slightly above Uq + ^ (see also the right panel of Fig.js}. Other con- 
tributions to the Fourier signal at the surface stem from the two edge 
QPOs e[ + ^ and E^ at frequencies just below Uq~^ and close 
to respectively. However, the edge QPOs are damped much 

faster than turning-point QPOs, and the amplitude in the Fourier 
analysis for the same initial data is about two orders of magni- 
tude smaller than for a reference model without crust. Therefore 
we doubt that edge QPOs are the explanation for the observed fre- 
quencies in SGR for magnetic fields B < 5 x 10 14 G. 

In Table [6] we give the maximum of the Fourier amplitude at 
the surface for a simulation with crust at 10 15 G with respect to the 
corresponding simulation without crust. The closer this value is to 
1 .0 the stronger is the amplitude of the QPO at the surface. For very 
low values the crust shields the QPOs efficiently, i.e. they are con- 
fined to the core of the neutron star. For the DH crustal EoS there is 
already a considerable amount of oscillations penetrating the crust 
and reaching the surface. We therefore argue that magneto-elastic 
oscillations break through the crust around 10 15 G. However, for 
the NV EoS the amplitudes for models with crust never reach 10% 
of the values for models without crust, as the crust is more ex- 
tended in this case, and QPOs can break through the crust only 
for even stronger magnetic fields. Nevertheless, the threshold of 



Model 


^ohJIO 15 Gl 


#28H 7 ri0 15 Gl 


APR+DH 1.4 


5.0 


4.6 


APR+DH 1.8 


7.1 


6.6 


APR+DH 2.2 


11.0 


10.3 


L+DH 1.4 


4.1 


3.8 




0.0 


C 1 

J.L 


L+DH 2.2 


6.4 


6.7 


APR+NV 1.4 


5.1 


4.8 


APR+NV 1.8 


7.4 


6.9 


APR+NV 2.2 


11.3 


10.6 


L+NV 1.4 


4.6 


4.4 


L+NV 1.8 


6.2 


5.8 


L+NV 2.2 


8.1 


7.6 



Table 7. Mean surface magnetic field strength required to match the fre- 
quency of the fundamental Alfven QPO to 30 Hz observed in SGR 1 806-20 
and 28 Hz of SGR 1900+14. 

B rsj 10 15 G should be a good approximation for all models. In- 
terestingly, this result is comparable with estimates of the magnetic 
field strengths for SGRs showing giant flares. 



6 SUMMARY AND DISCUSSION 

We have presented results from 2-dimensional, general-relativistic, 
magneto-hydrodynamical simulations of neutron stars with an ex- 
tended solid crust, in continuation of our initial results communi- 
cated as a letter ( Ga bler et al.|2011| K Performing a comprehensive 
set of simulations for several neutron star models and EoS, we have 
been able to confirm our previous findings regarding the QPO struc- 
ture for three different regimes of the magnetic field strength. Our 
main results can be summarized as follows: 

• For weak magnetic fields, B < 5 x 10 13 G, purely shear os- 
cillations of the crust dominate the evolution in the latter. For in- 
termediate magnetic field strengths, 5 x 10 13 < B < 10 15 G, the 
n = crustal modes are damped very efficiently into the core of 
the neutron star on timescales of a fraction (~ 0.04) of the Alfven 
crossing time of the star. For example, a model with B = 10 14 G 
has a damping timescale of r < 100 ms. This effectively rules out 
purely shear oscillations of the crust as a possible explanation for 
the observed QPOs in SGRs for the poloidal magnetic field config- 
urations studied here, since the observed QPOs survive for tens of 
seconds at estimated magnetic field strengths B > 6 x 10 14 G. We 
find a damping timescale dependence on the magnetic field which 
scales as ~ dominated by the ability of the Alfven continuum 
to absorb the energy of crustal modes. 

• In comparison to the n = modes, the n > modes of the 
crust have damping timescales of order hundreds of milliseconds 
at around 5 x 10 14 G. The spatial structure of the n > 1 modes 
becomes significantly distorted in the presence of such a strong 
magnetic field. However, predictions for even stronger magnetic 
fields, B > 5 x 10 14 G, are currently not possible, because the grid 
resolution needed to couple the n > shear modes to Alfven oscil- 
lations in the core is too high to perform simulations in reasonable 
time. We thus cannot safely exclude the n > crustal shear modes 
as possible explanation for the high frequency QPOs observed in 
SGR 1806-20 at 625 and 1840 Hz, if the magnetic field is a purely 
dipolar one. 

• For magnetic fields, B > 5 x 10 14 G, we find no sign of the 
existence of discrete crustal shear modes. This is in contrast to the 
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results of the model of van Hoven & Levin (2011), who proposed 
that weakly damped global QPOs, resembling purely crustal shear 
modes in the zero magnetic field limit, may survive if the corre- 
sponding shear mode frequency lies in-between adjacent continua. 
In the most promising of our models we find gaps between the first 
four continua, but even when choosing the equilibrium model so 
that the fundamental crustal mode frequency falls in one of the ex- 
isting gaps, we found that the mode was damped very efficiently 
by the coupling to the Alfven continuum in the core. Furthermore, 
at dipolar magnetic field strengths B > 10 15 G, necessary to fit 
the frequencies of crustal shear modes into a gap, we expect that 
the magnetic field begins to have significant influence on the shear 
oscillations of the crust ( [Messios et al.||200l] |Sotani et al.| [ 008a 
Shaisulta nov & Eich ler 2009). Our findings strongly support the 
interpretation that for high magnetic field strengths magnetar os- 
cillations are a strongly-coupled magneto-elastic system, where a 
division into purely crustal modes and Alfven oscillations is no 
longer valid. 

• In the intermediate magnetic field regime (5 x 10 13 G< B < 
10 15 G) the QPOs are largely confined to the core of the neutron 
star. We find three families of QPOs: upper, edge and lower QPOs. 
Their spatial structure coincides very well with the expectations 
from our semi-analytic model, if we assume that the oscillations are 
reflected at the crust-core boundary. Together with the strong damp- 
ing of the crustal shear modes, the reflection of Alfven QPOs at the 
crust-core interface leads to very small oscillation amplitudes in the 
crust. Moreover, when changing the magnetic field strength, the po- 
sition of the maximum of the corresponding upper QPO within the 
star changes significantly due to the interaction with the crust. 

• We have also determined the dipolar magnetic field strength at 
which the magneto-elastic QPOs break through the crust and reach 
the surface with significant amplitudes. This happens around B ~ 
10 15 G for the DH crust EoS and at slightly stronger magnetic fields 
for the NV crust EoS. This difference can be understood easily, 
because the NV EoS leads to thicker crusts and larger shear moduli, 
in particular near the crust-core interface. 

• Between 10 15 G< B < 5 x 10 15 G the QPOs are still trans- 
forming from being confined to the core to being able to reach the 
surface and the QPOs have spatial structures different from those 
at much lower or much higher magnetic field strengths. Colaiuda 
|& Ko kkotas (2011) report global, discrete Alfven modes in gaps 
between continua, at magnetic fields around 4 x 10 15 G. In our 
model, this is still in the above transition regime, and the reported 
oscillations are not discrete Alfven modes, but rather an effect of 
the transition between the two limiting regimes. 

• For dipolar magnetic field strengths B > 5 x 10 15 Gthe 
magneto-elastic oscillations have an almost Alfven-like character 
in the whole star and the role of the shear modulus in the crust is 
diminished, recovering the results of Sotani et al. (008b) and Cerda- 
|Duran et al.| ( [2009l >. 

The model we have presented allows for a tentative interpre- 
tation of the observed QPOs in terms of the predominantly Alfven 
QPOs which reach the surface. The main family of QPOs we expect 
to play a role in the observed QPOs are the upper (turning-point) 
QPOs. As first pointed out by |Sotani et al.| ( |008a| ), the interpre- 
tation of the frequencies 30, 92 and 150 Hz in SGR 1806-20 in 
terms of the Uq~\ and QPOs is very tempting. Here 
we note that the 18 Hz oscillation may be interpreted as the first 
edge QPO Eq~\ which for model APR+DH has a frequency of 

about 0.57 x /,,(-) ~ 17 Hz (this could also correspond to the fre- 

■ ■ i 

quency at 16.9 Hz found by Hamba ryan et al.|2011| ). Similarly, the 



36.4 Hz oscillation could be interpreted as the E{~ } QPO. How- 
ever, one should be cautious with these latter identifications, since 
edge QPOs are not as long-lived as turning-point QPOs. 

If we require that the fundamental upper QPO matches the 
observed 30 Hz QPO in SGR 1806-20 or the 28 Hz QPO in SGR 
1900+14, we obtain mean dipolar magnetic field strengths at the 
surface in the possible range 3.8 x 10 15 G< B < 1.1 x 10 16 G 
(for the particular choices of EoS and masses) as reported in Table 
[7] This is a rather narrow range and it is at only somewhat larger 
magnetic field strengths than simple estimates for magnetic fields 
in known magnetar s. 

We conclude by highlighting that our model provides two dif- 
ferent constraints on the magnetic field strength. The first con- 
straint is that the dipolar magnetic field strength must be larger 
than B ^ 10 15 G for QPOs to break out of the crust, which is a 
lower limit on the magnetic field. This constraint is independent 
of a particular identification of observed QPOs and only depends 
on the assumed magnetic field structure of a pure dipole. The sec- 
ond constraint comes from the matching of the lowest-frequency 
observed QPO, that appears at near-integer multiples, with the fun- 
damental Uq~^ QPO. This constrains the mean dipolar magnetic 
field strength at the surface of the neutron star to be in the range 
of 3.8 x 10 15 G< B < 1.1 x 10 16 G for the sample of EoSs and 
various masses that we assumed here. 

The above constraints, favoring somewhat stronger magnetic 
fields than estimated for known magnetars, hint at a possible de- 
viation of the actual magnetic field structure from a global dipole. 
This is not surprising, as it is known that a purely dipolar or purely 
toroidal field is not a stable magnetic field configuration in com- 
pact stars (see |Braithwaite & Nordlundl|2006| |Lasky et al.]|2011| 
ICiolfi et al.|2011||Kiuchi et al.|2011| and references therein). We 
may thus have an observational indication that the structure of the 
magnetic field in magnetars is in fact more complicated than a pure 
dipole. Other physical effects which may cause the magnetic field 
required to match the magneto-elastic QPOs to the range of ob- 
served frequencies to be weaker, are superfluidity of the neutrons 
and superconductivity of the protons in the core. 

We are planning to investigate the effects produced by chang- 
ing the magnetic field configuration and those caused by superflu- 
idity and superconductivity of the neutrons and protons in the core, 
which are expected to influence the Alfven speeds and the over- 
all dynamics (Passamonti & Andersson 2012). For the former, we 
need to investigate different dipolar configurations, mixed poloidal- 
toroidal configurations, relax the assumption of axisymmetry and 
include the coupling to polar oscillations. Moreover, the coupling 
to an exterior magnetosphere, where the X-ray emission is modu- 
lated, also has to be included in a complete model. 
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with u v and obtain 



(A6) 



APPENDIX A: ALTERNATIVE NUMERICAL METHOD 

In this appendix we derive a wave equation to describe the cou- 
pled crust-core oscillations of magnetars. To this end we expand 
all dynamical variables /(r, t) into a static unperturbed part /(r), 
denoted by a caret, and a time-dependent perturbation £/(r, £): 



f(T,t) = f(r) + 5f(T,t). 



(Al) 



For clarity we will omit the arguments r and t from now on. 

In the subsequent subsections we derive the equation for the 
displacement, describe the numerical implementation and compare 
the performance of this alternative method with the standard ap- 
proach using Riemann solvers on both sides of the crust-core inter- 
face. 



(Al) 



where we have used u v u v — — 1 and u u b u = . From u v b v ^ 



(see Papadopoulos & Esposito 1982) it follows that 



(A8) 



When linearizing Eq. ( |A6| > and using ( |A8| > Faraday's equation be- 
comes: 



u^Stfp = —b u .^5u u + h uf3 (5uf3-\b x + up-\5b x ) 
-\-u v Su^b x up-\ — 5u^J) v — u^5b u 



-\-u v bP (5u@-\u x + Uf3-\5u x ) 
-\-up-\u x (b^ ] 8u u + Sb^u") . 



(A9) 



Taking the (^-component of this equation, we arrive at a relation be- 
tween Sb^ and the spatial derivatives of the displacement £^ which 
reads 



Al Linear wave equation 

To derive the linearized wave equation governing the evolution 
of the displacement in the crust of a magnetized neutron star, we 
project the conservation equation of energy-momentum 



h%T% 



(ph + b 2 )u^ u u v 



0, 

+h»Jb?b» 



(A2) 



(A3) 



with h^ u = g^ v ' +u^u v and apply the following simplifications: (i) 
we linearize in the perturbations Sf and (ii) neglect all metric per- 
turbations (Cowling approximation, 5g = 0). Then Eq. \A2\ reads 

(ph + b 2 )Su^ u u u — 

— (dph + p5h + 2bp5b^ u^ v u — (ph + b 2 )u^ v 5u v 



b p b v 



■sT 



+h^ [b?5b u + SbPb" - f v (5P + bpSb?)] 



-2K 



(A4) 



Next we restrict ourselves to (iii) axisymmetry (f j(p = 0), and (iv) 
axial perturbations (8u l Su r 5u° 0, 5b r Sb e 0, 
and Sps — dh — Sp = SP — 0). Furthermore, we consider (v) 
the spherical symmetric, non-rotating background described by the 
line element ds 2 = -a 2 dt + <t 4 (dr 2 + r 2 d0 2 + r 2 sin(0) 2 d<p 2 ), 
and (vi) purely poloidal background fields (b* = 0). Applying all 
these simplifications one arrives at the following equation for Su^: 



(ph + b )u t 5u p t 



+5b* 



2$ 



b r 6b* r + b 9 6b* e - 2p s S^ (s) 

r - + - + ) V + 2cot(0)&' 
r a 1 



(A5) 



Because of the dependence of 1T V on the displacement we ex- 
press all other quantities in terms of the latter. Recalling the defini- 
tion of the corresponding velocity £^ = Sv? /u 1 ; , see Eq. (24i, the 
perturbed magnetic field Sb^ remains the only missing ingredient. 

To find an expression relating Sb^ to we contract the Fara- 
day equation 



\ — — r-fr r Su^ + b r Su^r + b Su^n 



+b r u\ r fj? t , 



where we have used that u* = a -1 . Plugging this relation into 
Eq. \A5\ we obtain 



(A10) 
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with 

Ao = (ph + b 2 ) (u* f , 
A 1 = A 3 b r + b r b\ + b e b r e + 1 Ms 
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and 

A 3 



4 ?V + ? + jr + 2 cot( ^)^ . (A16) 



The solution of the coefficient determinant of the second-order 
derivatives of Eq. \A\2\ 
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leads to the following eigenvalues: 



Ai 
A 2 

A 3 



1, 

As + b 2 
a 2 A 

As 



& 2 Ao 



(A19) 
(A20) 

(A21) 



where the eigenvector corresponding to A2 (A3) is oriented along 
(perpendicular to) the magnetic field lines, i.e. they correspond 
to the linearized version of Eq. ( [37] ). All eigenvalues are real, and 
hence Eq. \A\2\ is hyperbolic. 
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A2 Numerical implementation of the linearized wave 
equation in the crust 

As when using Riemann solvers on both sides of the crust-core 
interface, the equations and the numerical scheme in the fluid core 
of the neutron star remain unmodified. In order to evolve Eq. \A\2\ 
numerically in the crust we split it into two equations for £^ and 

and then perform an explicit Runge-Kutta integration. Since 
we are now evolving two systems with different variables in the 
core and in the crust, we have to impose interface conditions. 

The evolutions computed with this method rapidly become un- 
stable when increasing the magnetic field strength. It was therefore 
necessary to add some artificial dissipation. We used a fourth-order 
Kreiss-Oliger term eDT>4f, where V^f is the fourth-order numer- 
ical derivative of any function /. The minimal coefficient found to 
give stable evolutions is cd = 10 -2 . We checked the code to en- 
sure that this additional term does not influence the results of the 
simulations significantly. 



A3 Comparison of the two methods 

To compare the results of the crustal mode damping (see Sec- 
tion |5.2|) o btained with the two numerical methods presented in 
Section[2] [4] and in Appendix |A| we plot the evolution of the ve- 
locity at some point in the crust near the pole (Fig. |Al| ). Without 
magnetic field, the linear method is less dissipative, which is proba- 
bly related to the set-up of the interface conditions at the crust-core 
interface in this particular case. While we use the general condi- 
tions described in Section]!] for the Riemann solver approach, it 
is possible to use a simplified expression for the linear method. 
Because there is no magnetic field the expression of continuous 
traction at the crust-core interface leads to £f r = as at the sur- 
face. This provides a source of dissipation for the Riemann solver 
method, but none for the linear approach. In the presence of mag- 
netic fields the picture changes, and the evolution computed with 
the Riemann solver has less dissipation of crustal modes than the 
linear method (see Fig. |Al| ). When using different coefficients in 
the Kreiss-Oliger term, the curves for the linear method are indis- 
tinguishable. We can therefore rule out that the artificial dissipation 
dominates the numerical damping observed for the linear approach. 

Because of its superior behavior in the more generic case we 
generally used the Riemann solver method to obtain the numerical 
results. We checked that the linear approach agrees on the extracted 
frequencies with the Riemann solver for all regimes of the magnetic 
field strength considered by us. 



Figure Al. Evolution of the velocity f ,t at some point in the crust near the 
pole for B = 5 x 10 13 G and without a magnetic field. The two different 
numerical methods are denoted by lin for the linear method and R for the 
Riemann solver approach. The linear method is less dissipative for zero 
magnetic field, while the opposite holds when the magnetic field is turned 
on. 



APPENDIX B: EIGENMODES IN THE CRUST 



Without magnetic field Eq. \A\2\ simplifies to 



with 
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Assuming a harmonic time dependence ^(r, t) = ^(r)e %U}t and 
a separation of variables ^(r, 0) = R(r)G(0) leads to 
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where a prime denotes the derivative with respect to the corre- 
sponding variable r or 0. The angular and radial part have to fulfill 
this equation independently, such that: 

= 6 // (6') + 3cot6'e / (6') + A0e(6'), 



= ^i?"(r) + 
$4 



$4 



(B7) 
+ X 2 r R(r),(BS) 



where = (uj 2 Ao + . Both equations are of singular Sturm- 
Liouville type, and therefore their solutions R\ r and <d\ 9 form a 
complete, orthonormal set £;(r, 0) — R\ r G\ d . The solution of 
the angular part consists of the angular part of the vector spherical 
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harmonics ty n (6), which is related to the Legendre polynomials 

Pn(0)by 



dP n (0) 

80 



(B9) 



The equation to obtain R(r) is solved numerically with a shooting 
method. Because the eigenfunctions Hi (r, 6) form a complete set, 
it is possible to expand any displacement in terms of the former: 



(BIO) 



where the eigenmode coefficients Ai(t) are given by the inner 
product 



Mt) = (e(r,9,t),E t (r,d)) 



f 



^(r,e,t) Zi(r,6)wew r drd6 . 



(Bll) 



The corresponding weighting functions iu r and we are, according 
to the Sturm-Liouville theory 



we — sin(0) , 



(B12) 
(B13) 



aCrcc)- 1 ^^) 10 r4 cP (r cc )/i(r cc ) " 

We calculate the overlap integrals defined in Eq. \B\\\ with a 
fourth-order Simpsons rule algorithm. 
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